1 Setup

Code

library(easyRasch) # devtools::install_github("pgmj/easyRasch")
#library(RISEkbmRasch)
library(grateful)
library(ggrepel)
library(car)
library(kableExtra)
library(tidyverse)
library(eRm)
library(iarm)
library(mirt)
library(psych)
library(psychotree)
library(matrixStats)
library(reshape)
library(knitr)
library(patchwork)
library(formattable) 
library(glue)


### optional libraries
#library(TAM)
#library(skimr)
#library(janitor)
# 935600
# nohup quarto render MADRS.qmd --to html &> render_madrs.out &

### some commands exist in multiple packages, here we define preferred ones that are frequently used
select <- dplyr::select
count <- dplyr::count
recode <- car::recode
rename <- dplyr::rename

path_prim <- '/srv/projects/nils/study4rasch'
path2 <- '/Volumes/projects/k8_CPF_Kaldo/Data/Lärande Maskiner/Nils/Study4data'
running_machine = Sys.info()[['sysname']]
path = ifelse(running_machine!='Linux',path2,path_prim)
set_mywd_path = ifelse(running_machine!='Linux','./study4rasch/rasch_measures','./rasch_measures')
n_cores =ifelse(running_machine!='Linux',1,8)

tryCatch({setwd(set_mywd_path)}, error = function(set_mywd_path){print('already in right directory')})

[1] "already in right directory"

Code

source('settings.R') # containing item labels 
source('mod_easyRasch_func.R') # modified functions + added

2 Importing data

Code

### import data
df_raw <- read.csv(file=file.path(path,'data','raw_data','outcome.csv')) #

### Clean data 
source("preprocess.R")
df <- preprocess_data(df_raw)

# items to use 
items_to_use <- 'MADRS.S' # options are #MADRS.S, LSAS.SR, PDSS.SR
itemlabels <- itemlabels %>% 
  as_tibble() %>%
  filter(str_detect(itemnr, items_to_use))


# DIF+aux variables to use 
# DIF not used 'marital_status','children', 'working'
dif_variables <- c('Patient','Treatment','sex','age','TreatmentAccessStart','education')

### Make a backup of the dataframe, in case you need to revert changes at some point
df.all <- df

print(itemlabels)

# A tibble: 9 × 2
  itemnr    item                      
  <chr>     <chr>                     
1 MADRS.S_1 Sadness                   
2 MADRS.S_2 Anxiety                   
3 MADRS.S_3 Reduced sleep             
4 MADRS.S_4 Reduced appetite          
5 MADRS.S_5 Concentration difficulties
6 MADRS.S_6 Lassitude                 
7 MADRS.S_7 Inability to feel         
8 MADRS.S_8 Pessimism                 
9 MADRS.S_9 Suicidal thoughts

3 Checking missing

Code

itemStart <- items_to_use
out = RImissing_mod(df,itemStart,facet="Time")
print(out)

3.1 Remove missing

Code

##### Before filtering out participants, you should check the missing data structure using RImissing() and RImissingP()

# If you want to include participants with missing data, input the minimum number of items responses that a participant should have to be included in the analysis:
min.responses <- 4 # In our case if they have missing on one item all items for that time are missing. 

# Select the variables we will work with, and filter out respondents with missing data
df <- df %>% 
  select(all_of(c(dif_variables,"Time",itemlabels$itemnr))) %>%  # v
  filter(rowSums(is.na(select(.,all_of(itemlabels$itemnr)))) < min.responses)

4 Create DIF variables

Code

#---- Create DIF variables----
  
# DIF variables into vectors, recoded as factors since DIF functions need this
# these could also be stored in its own dataframe (not a tibble) instead of as vectors
# Named vector for the new types

type_transform <- c(Treatment = "factor", sex = "factor", age = "numeric",
                    TreatmentAccessStart ="integer",education="factor",Time='factor')

# Transform columns based on the named vector
dif <- df %>%
  mutate(across(names(type_transform), ~ switch(type_transform[which(names(type_transform) == cur_column())], 
                                               "integer" = as.integer(.), 
                                               "character" = as.character(.),
                                               "factor" = as.factor(.), 
                                               "numeric" = as.numeric(.)))) %>%
                                               as.data.frame(.) %>%
                                               select(!all_of(itemlabels$itemnr))


# then remove them from dataframe, since we need a dataframe with only item data for the Rasch analyses
df <- df %>% select(all_of(c("Time",itemlabels$itemnr))) # add time here ? 
dfnotime <- df %>% select(!Time)
source("RISE_theme.R")

5 Demographics descriptives

Code

dif_spec <- dif %>% filter(Time=='SCREEN') %>% select(!all_of(c("Time","Patient")))
summary(dif_spec)
# RIdemographics(dif_spec,'Treatment') <- function crashes TODO make own

 Treatment  sex           age        TreatmentAccessStart         education   
 MDD:2987   F:4005   Min.   :18.00   Min.   :2008         Primary      : 504  
 PD :1719   M:2456   1st Qu.:27.00   1st Qu.:2012         Secondary    :3152  
 SAD:1755            Median :33.00   Median :2014         Postsecondary:2805  
                     Mean   :35.38   Mean   :2014                             
                     3rd Qu.:42.00   3rd Qu.:2017                             
                     Max.   :84.00   Max.   :2019

6 Overall number of responses

Code

# Collapsed 
RIallresp(dfnotime)
# Seperate 
RIallresp_over_times <- df %>% # 
  split(.$Time) %>% # split the data 
  map(~ RIallresp(.x %>% dplyr::select(!Time))) #+ labs(title = .x$Time)) # create separate for each time

# make nice later 
RI_allresp_kable_grid = combine_kables_grid(RIallresp_over_times,cols=3)
RI_allresp_kable_grid

Response category	Number of responses	Percent
0	147379	24.9
1	148593	25.1
2	152831	25.8
3	67964	11.5
4	61598	10.4
5	10877	1.8
6	2301	0.4

SCREEN.Response category	SCREEN.Number of responses	SCREEN.Percent	PRE.Response category	PRE.Number of responses	PRE.Percent	WEEK01.Response category	WEEK01.Number of responses	WEEK01.Percent
0	8232	14.2	0	9700	17.3	0	9626	20.1
1	7983	13.7	1	9965	17.8	1	10868	22.7
2	15050	25.9	2	15616	27.9	2	13754	28.8
3	10126	17.4	3	9146	16.3	3	6674	14.0
4	13231	22.8	4	9629	17.2	4	5814	12.2
5	2816	4.8	5	1608	2.9	5	935	2.0
6	711	1.2	6	307	0.5	6	155	0.3

WEEK02.Response category	WEEK02.Number of responses	WEEK02.Percent	WEEK03.Response category	WEEK03.Number of responses	WEEK03.Percent	WEEK04.Response category	WEEK04.Number of responses	WEEK04.Percent
0	10426	21.2	0	10968	22.7	0	11626	24.7
1	11747	23.9	1	12365	25.6	1	12516	26.6
2	14076	28.6	2	13422	27.8	2	12690	26.9
3	6367	12.9	3	5755	11.9	3	5254	11.2
4	5502	11.2	4	4697	9.7	4	4155	8.8
5	884	1.8	5	856	1.8	5	722	1.5
6	183	0.4	6	150	0.3	6	134	0.3

WEEK05.Response category	WEEK05.Number of responses	WEEK05.Percent	WEEK06.Response category	WEEK06.Number of responses	WEEK06.Percent	WEEK07.Response category	WEEK07.Number of responses	WEEK07.Percent
0	11754	26.2	0	11745	27.3	0	11877	28.9
1	12299	27.4	1	12301	28.6	1	11965	29.1
2	11657	26.0	2	10933	25.4	2	10044	24.5
3	4732	10.6	3	4201	9.8	3	3801	9.3
4	3639	8.1	4	3242	7.5	4	2804	6.8
5	639	1.4	5	525	1.2	5	460	1.1
6	100	0.2	6	136	0.3	6	98	0.2

WEEK08.Response category	WEEK08.Number of responses	WEEK08.Percent	WEEK09.Response category	WEEK09.Number of responses	WEEK09.Percent	WEEK10.Response category	WEEK10.Number of responses	WEEK10.Percent
0	12015	30.6	0	11977	32.1	0	11951	34.0
1	11516	29.3	1	11197	30.0	1	10795	30.7
2	9632	24.5	2	8719	23.4	2	7755	22.0
3	3218	8.2	3	2913	7.8	3	2580	7.3
4	2431	6.2	4	2113	5.7	4	1775	5.0
5	408	1.0	5	335	0.9	5	272	0.8
6	92	0.2	6	69	0.2	6	71	0.2

POST.Response category	POST.Number of responses	POST.Percent
0	15482	34.9
1	13076	29.5
2	9483	21.4
3	3197	7.2
4	2566	5.8
5	417	0.9
6	95	0.2

7 Descriptives of raw data

Response distribution for all items are summarized below.

7.1 Floor/ceiling effects

Code

RIrawdist(dfnotime)

Code

#df_test = df %>% filter(Time=='SCREEN') %>% select(!Time)
# Across time 

# seperated time - might need to fix this remake in ggplot?
png(filename=paste0("./plots/",itemStart,"_rawdist_over_time.png"),
width=800,height=2400)
par(mfrow=c(7,2),mar = c(4, 4, 4, 2))  # c(5 <- bottom lines, 4, 4, 2)
RIrawdist_over_times <- df %>% # 
  split(.$Time) %>% # split the data 
  imap(~ {
    RIrawdist(.x %>% dplyr::select(!Time))
    mtext(paste("Time:", .y), side = 1, line = 3, adj = 0, cex = 0.8)
    })

dev.off()

png 
  2

Code

print('Saved plot to /plots')

[1] "Saved plot to /plots"

7.2 Guttman structure

Code

RIheatmap(dfnotime)

Code

ggplots_with_spec_func_over_time(df,func_call = RIheatmap,title_all='Guttman structure over time',blankx=TRUE)

7.3 Descriptives - item level

Code

# over all times 
RItileplot(dfnotime)

Code

# by time 
ggplots_with_spec_func_over_time(df,func_call = RItileplot,title_all='',blankx=FALSE)

Code

RIbarstack(dfnotime)

Code

# by time 
ggplots_with_spec_func_over_time(df,func_call = RIbarstack,title_all='',blankx=FALSE)

Code

RIbarplot(dfnotime)

8 Rasch analysis 1

The eRm package, which uses Conditional Maximum Likelihood (CML) estimation, will be used primarily. For this analysis, the Partial Credit Model will be used.

Code

df_test = df %>% filter(Time=='SCREEN') %>% select(!Time)
RIlistItemsMargin(df, fontsize = 12)

itemnr	item
MADRS.S_1	Sadness
MADRS.S_2	Anxiety
MADRS.S_3	Reduced sleep
MADRS.S_4	Reduced appetite
MADRS.S_5	Concentration difficulties
MADRS.S_6	Lassitude
MADRS.S_7	Inability to feel
MADRS.S_8	Pessimism
MADRS.S_9	Suicidal thoughts

8.1 Timepoint decision

Due to the amount of data across timepoints targeting will first be inspected, and the best fitting timepoint will be chosen based on targeting, which in turn will inform the rest of the analyses.

8.1.1 Targeting

Code

# increase fig-height above as needed, if you have many items og. value was 5. 
#RItargeting(dfnotime)

require(patchwork)
targeting_time_plots <- ggplots_with_spec_func_over_time(
  df,func_call = RItargeting,title_all='',
  blankx=FALSE,colsplot = 3,rowsplot=5,return_list_plots=TRUE)
good_layout_targeting_plots <- wrap_plots(targeting_time_plots, ncol = 3) + plot_layout(guides= "collect")
#print(good_layout_targeting_plots)
ggsave(
  filename = paste0("./plots/",items_to_use,"_targeting_over_time.png"),
  plot = good_layout_targeting_plots,
  width = 20,  # Increase width (in inches)
  height = 45, # Increase height (in inches)
  dpi = 300    # High resolution
)
print(good_layout_targeting_plots)

8.1.2 DIF treatment group targeting

Week02 looks promising.

Looking at treatment group for timepoints

Code

screenplot <- treatment_dif_by_time(timepoint='SCREEN',df=df,dif=dif,given_func=RItargeting)
wrap_plots(screenplot, ncol = 3)

Code

preplot <- treatment_dif_by_time(timepoint='PRE',df=df,dif=dif,given_func=RItargeting)
wrap_plots(preplot, ncol = 3)

Code

week01plot <- treatment_dif_by_time(timepoint='WEEK01',df=df,dif=dif,given_func=RItargeting)
wrap_plots(week01plot, ncol = 3)

Code

week02plot <- treatment_dif_by_time(timepoint='WEEK02',df=df,dif=dif,given_func=RItargeting)
wrap_plots(week02plot, ncol = 3)

8.2 Analysis continuation week02

Code

df_week02 <- df %>% filter(Time=='WEEK02') %>% select(!Time) #%>% sample_n(800)

Code

RIitemfit(df_week02, cutoff = "Smith98")

Item	InfitMSQ	Location	1 ± 2 / √n
MADRS.S_1	0.688	0.26	[0.973, 1.027]
MADRS.S_2	1.108	-0.78	[0.973, 1.027]
MADRS.S_3	1.593	-0.16	[0.973, 1.027]
MADRS.S_4	1.322	0.99	[0.973, 1.027]
MADRS.S_5	0.938	-0.17	[0.973, 1.027]
MADRS.S_6	0.793	-0.27	[0.973, 1.027]
MADRS.S_7	0.735	0.14	[0.973, 1.027]
MADRS.S_8	1.058	-0.53	[0.973, 1.027]
MADRS.S_9	0.846	0.60	[0.973, 1.027]
Note:
MSQ values based on conditional estimation (n = 5465 complete cases).

Code

simfit1 <- RIgetfit(df_week02, iterations = 300, cpu = n_cores)  

RIitemfit(df_week02, simfit1)

Item	InfitMSQ	Infit thresholds	OutfitMSQ	Outfit thresholds	Infit diff	Outfit diff	Location
MADRS.S_1	0.688	[0.95, 1.045]	0.677	[0.953, 1.054]	0.262	0.276	0.26
MADRS.S_2	1.108	[0.954, 1.046]	1.118	[0.953, 1.05]	0.062	0.068	-0.78
MADRS.S_3	1.593	[0.955, 1.052]	1.71	[0.939, 1.06]	0.541	0.650	-0.16
MADRS.S_4	1.322	[0.951, 1.056]	1.437	[0.923, 1.113]	0.266	0.324	0.99
MADRS.S_5	0.938	[0.955, 1.045]	0.929	[0.957, 1.051]	0.017	0.028	-0.17
MADRS.S_6	0.793	[0.954, 1.041]	0.788	[0.954, 1.042]	0.161	0.166	-0.27
MADRS.S_7	0.735	[0.955, 1.04]	0.725	[0.955, 1.04]	0.220	0.230	0.14
MADRS.S_8	1.058	[0.95, 1.052]	1.083	[0.946, 1.056]	0.006	0.027	-0.53
MADRS.S_9	0.846	[0.965, 1.049]	0.803	[0.955, 1.062]	0.119	0.152	0.60
Note:
MSQ values based on conditional calculations (n = 5465 complete cases). Simulation based thresholds from 300 simulated datasets.

Code

RIgetfitPlot(simfit1, df_week02)

Code

ICCplot(as.data.frame(df_week02), 
        itemnumber = 1:4, 
        method = "cut", 
        itemdescrip = itemlabels$itemnr[1:4])

[1] "Please press Zoom on the Plots window to see the plot"

Code

### also suggested:
# library(RASCHplot) # install first with `devtools::install_github("ERRTG/RASCHplot")`
# CICCplot(PCM(df),
#          which.item = c(1:4),
#          lower.groups = c(0,7,14,21,28,35),
#          grid.items = TRUE)

Code

ICCplot(as.data.frame(df_week02), 
        itemnumber = 5:8, 
        method = "cut", 
        itemdescrip = itemlabels$itemnr[5:8])

[1] "Please press Zoom on the Plots window to see the plot"

Code

RIrestscore(df_week02)

Item	Observed value	Model expected value	Absolute difference	Adjusted p-value (BH)	Statistical significance level	Location	Relative location
MADRS.S_1	0.75	0.61	0.14	0.000	***	0.26	1.66
MADRS.S_2	0.59	0.62	0.03	0.000	***	-0.78	0.62
MADRS.S_3	0.44	0.63	0.19	0.000	***	-0.16	1.24
MADRS.S_4	0.46	0.61	0.15	0.000	***	0.99	2.38
MADRS.S_5	0.65	0.62	0.03	0.000	***	-0.17	1.22
MADRS.S_6	0.72	0.62	0.10	0.000	***	-0.27	1.13
MADRS.S_7	0.74	0.60	0.14	0.000	***	0.14	1.54
MADRS.S_8	0.63	0.64	0.01	0.055	.	-0.53	0.87
MADRS.S_9	0.69	0.60	0.09	0.000	***	0.60	1.99

Code

RIbootRestscore(df_week02,cpu=n_cores,iterations = 500,samplesize=800) # samplesize=600

Item	Item-restscore result	% of iterations	Conditional MSQ infit	Relative average item location
MADRS.S_1	overfit	100.0	0.69	1.66
MADRS.S_2	no misfit	82.4	1.11	0.62
MADRS.S_2	underfit	17.6	1.11	0.62
MADRS.S_3	underfit	100.0	1.59	1.24
MADRS.S_4	underfit	100.0	1.32	2.38
MADRS.S_5	no misfit	57.8	0.94	1.22
MADRS.S_5	overfit	42.2	0.94	1.22
MADRS.S_6	overfit	100.0	0.79	1.13
MADRS.S_7	overfit	100.0	0.74	1.54
MADRS.S_8	no misfit	95.0	1.06	0.87
MADRS.S_9	overfit	98.8	0.85	1.99
Note:
Results based on 500 bootstrap iterations with a sample size of 800. Conditional mean-square infit based on complete responders only, n = 5465.

Code

RIpcmPCA(df_week02)

PCA of Rasch model residuals
Eigenvalues	Proportion of variance
1.70	22.5%
1.50	16.9%
1.31	14.5%
1.08	14.1%
0.97	10.4%

Code

simcor1 <- RIgetResidCor(df_week02, iterations = 500, cpu = n_cores)
RIresidcorr(df_week02, cutoff = simcor1$p99)

	MADRS.S_1	MADRS.S_2	MADRS.S_3	MADRS.S_4	MADRS.S_5	MADRS.S_6	MADRS.S_7	MADRS.S_8
MADRS.S_1
MADRS.S_2	-0.06
MADRS.S_3	-0.25	-0.12
MADRS.S_4	-0.16	-0.15	0
MADRS.S_5	-0.16	-0.08	-0.15	-0.13
MADRS.S_6	-0.05	-0.17	-0.26	-0.13	0.14
MADRS.S_7	0.15	-0.26	-0.24	-0.14	-0.06	0.09
MADRS.S_8	-0.09	-0.01	-0.26	-0.26	-0.2	-0.15	-0.09
MADRS.S_9	0.12	-0.21	-0.23	-0.12	-0.21	-0.11	0.09	0.07
Note:
Relative cut-off value is -0.046, which is 0.061 above the average correlation (-0.107). Correlations above the cut-off are highlighted in red text.

Code

RIpartgamLD(df_week02)

Item 1	Item 2	Partial gamma	SE	Lower CI	Upper CI	Adjusted p-value (BH)
MADRS.S_6	MADRS.S_5	0.349	0.019	0.310	0.387	0.000
MADRS.S_1	MADRS.S_7	0.338	0.021	0.296	0.380	0.000
MADRS.S_1	MADRS.S_9	0.334	0.021	0.292	0.375	0.000
MADRS.S_7	MADRS.S_1	0.325	0.021	0.283	0.367	0.000
MADRS.S_7	MADRS.S_9	0.315	0.022	0.271	0.358	0.000
MADRS.S_5	MADRS.S_6	0.305	0.020	0.266	0.344	0.000
MADRS.S_9	MADRS.S_1	0.294	0.022	0.252	0.337	0.000
MADRS.S_9	MADRS.S_7	0.294	0.022	0.250	0.337	0.000
MADRS.S_7	MADRS.S_6	0.291	0.021	0.249	0.333	0.000
MADRS.S_9	MADRS.S_8	0.277	0.020	0.239	0.316	0.000
MADRS.S_6	MADRS.S_7	0.272	0.022	0.230	0.314	0.000
MADRS.S_8	MADRS.S_9	0.214	0.020	0.174	0.253	0.000
MADRS.S_8	MADRS.S_2	0.172	0.019	0.135	0.209	0.000
MADRS.S_2	MADRS.S_8	0.162	0.018	0.126	0.198	0.000
MADRS.S_1	MADRS.S_8	0.13	0.020	0.089	0.170	0.000
MADRS.S_4	MADRS.S_3	0.119	0.020	0.081	0.158	0.000
MADRS.S_1	MADRS.S_2	0.11	0.021	0.068	0.152	0.000
MADRS.S_1	MADRS.S_6	0.104	0.022	0.061	0.147	0.000
MADRS.S_3	MADRS.S_4	0.104	0.019	0.066	0.143	0.000
MADRS.S_7	MADRS.S_5	0.095	0.022	0.052	0.139	0.001
MADRS.S_7	MADRS.S_8	0.085	0.021	0.043	0.126	0.004
MADRS.S_5	MADRS.S_3	0.072	0.019	0.035	0.109	0.009
MADRS.S_5	MADRS.S_2	0.07	0.020	0.031	0.110	0.031

Code

RIloadLoc(df_week02)

Code

mirt(df_week02, model=1, itemtype='Rasch', verbose = FALSE) %>% 
  plot(type="trace", as.table = TRUE, 
       theta_lim = c(-6,6))

Code

# for fewer items or a more magnified figure, use:
#RIitemCats(df)

Code

# increase fig-height above as needed, if you have many items
RItargeting(df_week02)

Code

RIitemHierarchy(df_week02)

Code

iarm::score_groups(as.data.frame(df_week02)) %>% 
  as.data.frame(nm = "score_group") %>% 
  dplyr::count(score_group)

  score_group    n
1           1 2752
2           2 2713

Code

dif_plots <- df_week02 %>% 
  add_column(dif = iarm::score_groups(.)) %>% 
  split(.$dif) %>% # split the data using the DIF variable
  map(~ RItileplot(.x %>% dplyr::select(!dif)) + labs(title = .x$dif))
dif_plots[[1]] + dif_plots[[2]]

Code

clr_tests(df_week02, model = "PCM")


Conditional Likelihood Ratio Tests:

        clr  df pvalue sig 
overall 1359 52 0       ***

Code

item_obsexp(PCM(df_week02))

Score group 1: 
          mean obs mean exp std.res sig
MADRS.S_1  0.719    0.828   -8.015  -- 
MADRS.S_2  1.828    1.792    2.338  +  
MADRS.S_3  1.214    1.018   12.700  ++ 
MADRS.S_4  0.339    0.241   10.175  ++ 
MADRS.S_5  1.186    1.191   -0.344     
MADRS.S_6  1.111    1.180   -4.877  -- 
MADRS.S_7  0.912    1.006   -6.777  -- 
MADRS.S_8  1.426    1.416    0.608     
MADRS.S_9  0.577    0.639   -4.778  -- 

Score group 2: 
          mean obs mean exp std.res sig
MADRS.S_1   2.430    2.323    5.836 ++ 
MADRS.S_2   3.293    3.329   -2.123 -  
MADRS.S_3   2.516    2.708  -10.544 -- 
MADRS.S_4   1.264    1.361   -5.422 -- 
MADRS.S_5   2.719    2.714    0.276    
MADRS.S_6   2.828    2.760    3.704 ++ 
MADRS.S_7   2.403    2.311    5.610 ++ 
MADRS.S_8   3.369    3.379   -0.559    
MADRS.S_9   1.897    1.835    4.007 ++

Code

grouping_based_on_score = score_groups(as.data.frame(df_week02))
partgam_DIF(dat.items = as.data.frame(df_week02),
            dat.exo = grouping_based_on_score)

       Item                     Var gamma  se pvalue padj.BH sig lower upper
1 MADRS.S_1 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
2 MADRS.S_2 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
3 MADRS.S_3 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
4 MADRS.S_4 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
5 MADRS.S_5 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
6 MADRS.S_6 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
7 MADRS.S_7 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
8 MADRS.S_8 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
9 MADRS.S_9 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA

Code

dif.sex <- dif %>% filter(Time=='WEEK02') %>% select(sex) %>% as.vector(.)
dif.sex <- dif.sex[[1]] #female = 1, male = 2 
RIdifTable(df_week02, dif.sex)

Item	2	3	Mean location	StDev	MaxDiff
MADRS.S_1	0.260	0.243	0.251	0.013	0.018
MADRS.S_2	-0.834	-0.723	-0.779	0.078	0.111
MADRS.S_3	-0.135	-0.254	-0.194	0.084	0.119
MADRS.S_4	0.907	1.155	1.031	0.175	0.248
MADRS.S_5	-0.155	-0.258	-0.207	0.073	0.103
MADRS.S_6	-0.196	-0.433	-0.314	0.168	0.237
MADRS.S_7	0.192	0.004	0.098	0.133	0.188
MADRS.S_8	-0.594	-0.439	-0.516	0.110	0.155
MADRS.S_9	0.554	0.705	0.629	0.106	0.150

Code

partgam_DIF(dat.items = as.data.frame(df_week02),
            dat.exo = dif.sex)

       Item     Var   gamma     se pvalue padj.BH  sig   lower   upper
1 MADRS.S_1 dif.sex  0.2073 0.0263 0.0000  0.0000  ***  0.1558  0.2588
2 MADRS.S_2 dif.sex -0.1617 0.0239 0.0000  0.0000  *** -0.2085 -0.1148
3 MADRS.S_3 dif.sex  0.0010 0.0228 0.9641  1.0000      -0.0438  0.0458
4 MADRS.S_4 dif.sex -0.1020 0.0274 0.0002  0.0018   ** -0.1558 -0.0482
5 MADRS.S_5 dif.sex  0.0548 0.0254 0.0312  0.2811       0.0049  0.1046
6 MADRS.S_6 dif.sex  0.0786 0.0260 0.0025  0.0225    *  0.0277  0.1296
7 MADRS.S_7 dif.sex  0.2151 0.0266 0.0000  0.0000  ***  0.1630  0.2672
8 MADRS.S_8 dif.sex -0.1738 0.0233 0.0000  0.0000  *** -0.2195 -0.1280
9 MADRS.S_9 dif.sex  0.0961 0.0276 0.0005  0.0045   **  0.0420  0.1502

Code

dif.age <- dif %>% filter(Time=='WEEK02') %>% select(age) %>% as.vector(.)
dif.age <- dif.age[[1]]
RIdifTable(df_week02, dif.age)

Item	3	4	6	8	9	Mean location	StDev	MaxDiff
MADRS.S_1	0.201	0.220	0.308	0.272	-0.077	0.185	0.152	0.385
MADRS.S_2	-0.873	-0.405	-0.902	-0.444	-0.773	-0.679	0.238	0.497
MADRS.S_3	0.039	0.050	-0.279	-0.356	-0.192	-0.148	0.185	0.406
MADRS.S_4	0.854	0.544	1.000	0.902	0.687	0.798	0.181	0.456
MADRS.S_5	-0.160	-0.210	-0.229	-0.204	0.135	-0.133	0.152	0.364
MADRS.S_6	-0.315	-0.240	-0.274	-0.285	-0.430	-0.309	0.073	0.190
MADRS.S_7	0.149	0.118	0.217	-0.003	0.262	0.149	0.102	0.265
MADRS.S_8	-0.484	-0.305	-0.564	-0.476	-0.167	-0.399	0.161	0.397
MADRS.S_9	0.589	0.228	0.722	0.595	0.555	0.538	0.185	0.494

Code

partgam_DIF(dat.items = as.data.frame(df_week02),
            dat.exo = dif.age)

       Item     Var   gamma     se pvalue padj.BH  sig   lower   upper
1 MADRS.S_1 dif.age -0.0380 0.0158 0.0161  0.1452      -0.0690 -0.0070
2 MADRS.S_2 dif.age -0.1415 0.0140 0.0000  0.0000  *** -0.1689 -0.1142
3 MADRS.S_3 dif.age  0.1160 0.0132 0.0000  0.0000  ***  0.0902  0.1418
4 MADRS.S_4 dif.age -0.0367 0.0158 0.0201  0.1808      -0.0676 -0.0058
5 MADRS.S_5 dif.age  0.0134 0.0149 0.3690  1.0000      -0.0158  0.0427
6 MADRS.S_6 dif.age -0.0317 0.0152 0.0369  0.3317      -0.0614 -0.0019
7 MADRS.S_7 dif.age  0.1092 0.0158 0.0000  0.0000  ***  0.0782  0.1401
8 MADRS.S_8 dif.age -0.0270 0.0137 0.0485  0.4366      -0.0539 -0.0002
9 MADRS.S_9 dif.age  0.0690 0.0161 0.0000  0.0002  ***  0.0374  0.1006

Code

dif.edu <- dif %>% filter(Time=='WEEK02') %>% select(education) %>% as.vector(.)
dif.edu <- dif.edu[[1]]
RIdifTable(df_week02, dif.edu)

Item	2	3	Mean location	StDev	MaxDiff
MADRS.S_1	0.271	0.210	0.241	0.043	0.061
MADRS.S_2	-0.792	-0.782	-0.787	0.007	0.010
MADRS.S_3	-0.159	-0.178	-0.169	0.014	0.020
MADRS.S_4	0.882	1.169	1.026	0.203	0.287
MADRS.S_5	-0.206	-0.154	-0.180	0.037	0.053
MADRS.S_6	-0.333	-0.190	-0.261	0.101	0.143
MADRS.S_7	0.144	0.109	0.126	0.025	0.035
MADRS.S_8	-0.454	-0.678	-0.566	0.158	0.224
MADRS.S_9	0.647	0.493	0.570	0.109	0.154

Code

partgam_DIF(dat.items = as.data.frame(df_week02),
            dat.exo = dif.edu)

       Item     Var   gamma     se pvalue padj.BH  sig   lower   upper
1 MADRS.S_1 dif.edu  0.0241 0.0238 0.3127  1.0000      -0.0226  0.0708
2 MADRS.S_2 dif.edu -0.0422 0.0215 0.0495  0.4453      -0.0844 -0.0001
3 MADRS.S_3 dif.edu  0.0231 0.0202 0.2529  1.0000      -0.0165  0.0627
4 MADRS.S_4 dif.edu -0.0631 0.0239 0.0082  0.0734    . -0.1099 -0.0164
5 MADRS.S_5 dif.edu -0.0262 0.0228 0.2487  1.0000      -0.0708  0.0183
6 MADRS.S_6 dif.edu -0.0888 0.0232 0.0001  0.0012   ** -0.1343 -0.0432
7 MADRS.S_7 dif.edu  0.0704 0.0243 0.0037  0.0333    *  0.0229  0.1180
8 MADRS.S_8 dif.edu  0.0871 0.0214 0.0000  0.0004  ***  0.0451  0.1290
9 MADRS.S_9 dif.edu  0.0530 0.0247 0.0320  0.2878       0.0046  0.1015

8.3 DIF-Year

Code

dif.yr <- dif %>% filter(Time=='WEEK02') %>% select(TreatmentAccessStart) %>% as.vector(.)
dif.yr <- dif.yr[[1]]
RIdifTable(df_week02, dif.yr)

Item	2	3	Mean location	StDev	MaxDiff
MADRS.S_1	0.329	0.213	0.271	0.082	0.116
MADRS.S_2	-0.778	-0.807	-0.792	0.021	0.029
MADRS.S_3	-0.172	-0.175	-0.174	0.002	0.003
MADRS.S_4	1.051	0.946	0.998	0.074	0.105
MADRS.S_5	-0.088	-0.240	-0.164	0.107	0.152
MADRS.S_6	-0.319	-0.253	-0.286	0.047	0.066
MADRS.S_7	0.066	0.182	0.124	0.083	0.117
MADRS.S_8	-0.557	-0.534	-0.546	0.016	0.023
MADRS.S_9	0.469	0.669	0.569	0.141	0.200

Code

partgam_DIF(dat.items = as.data.frame(df_week02),
            dat.exo = dif.yr)

       Item    Var   gamma     se pvalue padj.BH  sig   lower   upper
1 MADRS.S_1 dif.yr -0.0132 0.0165 0.4233  1.0000      -0.0455  0.0191
2 MADRS.S_2 dif.yr  0.0585 0.0151 0.0001  0.0010   **  0.0288  0.0881
3 MADRS.S_3 dif.yr -0.0031 0.0142 0.8297  1.0000      -0.0310  0.0249
4 MADRS.S_4 dif.yr -0.0271 0.0169 0.1084  0.9757      -0.0602  0.0060
5 MADRS.S_5 dif.yr  0.0094 0.0157 0.5488  1.0000      -0.0213  0.0402
6 MADRS.S_6 dif.yr -0.0044 0.0160 0.7812  1.0000      -0.0358  0.0269
7 MADRS.S_7 dif.yr -0.0177 0.0168 0.2921  1.0000      -0.0507  0.0152
8 MADRS.S_8 dif.yr  0.0228 0.0147 0.1222  1.0000      -0.0061  0.0517
9 MADRS.S_9 dif.yr -0.0704 0.0168 0.0000  0.0002  *** -0.1033 -0.0375

8.4 Analysis part 1 decision

Remove item 3 to begin with : motivation goes here

9 Rasch analysis 2

Code

df_r <- df_week02 %>% select(!MADRS.S_3)

Code

RIlistItemsMargin(df_r, fontsize = 12)

itemnr	item
MADRS.S_1	Sadness
MADRS.S_2	Anxiety
MADRS.S_4	Reduced appetite
MADRS.S_5	Concentration difficulties
MADRS.S_6	Lassitude
MADRS.S_7	Inability to feel
MADRS.S_8	Pessimism
MADRS.S_9	Suicidal thoughts

Code

RIitemfit(df_r, cutoff = "Smith98")

Item	InfitMSQ	Location	1 ± 2 / √n
MADRS.S_1	0.719	0.27	[0.973, 1.027]
MADRS.S_2	1.234	-0.88	[0.973, 1.027]
MADRS.S_4	1.501	1.06	[0.973, 1.027]
MADRS.S_5	1.03	-0.20	[0.973, 1.027]
MADRS.S_6	0.831	-0.30	[0.973, 1.027]
MADRS.S_7	0.771	0.14	[0.973, 1.027]
MADRS.S_8	1.131	-0.60	[0.973, 1.027]
MADRS.S_9	0.891	0.61	[0.973, 1.027]
Note:
MSQ values based on conditional estimation (n = 5465 complete cases).

Code

simfit1 <- RIgetfit(df_r, iterations = 300, cpu = n_cores)  

RIitemfit(df_r, simfit1)

Item	InfitMSQ	Infit thresholds	OutfitMSQ	Outfit thresholds	Infit diff	Outfit diff	Location
MADRS.S_1	0.719	[0.944, 1.048]	0.71	[0.94, 1.046]	0.225	0.23	0.27
MADRS.S_2	1.234	[0.954, 1.046]	1.254	[0.954, 1.049]	0.188	0.205	-0.88
MADRS.S_4	1.501	[0.943, 1.059]	1.811	[0.906, 1.179]	0.442	0.632	1.06
MADRS.S_5	1.03	[0.96, 1.039]	1.023	[0.959, 1.038]	no misfit	no misfit	-0.20
MADRS.S_6	0.831	[0.952, 1.048]	0.832	[0.953, 1.047]	0.121	0.121	-0.30
MADRS.S_7	0.771	[0.956, 1.046]	0.758	[0.953, 1.051]	0.185	0.195	0.14
MADRS.S_8	1.131	[0.951, 1.054]	1.16	[0.95, 1.056]	0.077	0.104	-0.60
MADRS.S_9	0.891	[0.951, 1.04]	0.842	[0.942, 1.054]	0.06	0.1	0.61
Note:
MSQ values based on conditional calculations (n = 5465 complete cases). Simulation based thresholds from 300 simulated datasets.

Code

RIgetfitPlot(simfit1, df_r)

Code

ICCplot(as.data.frame(df_r), 
        itemnumber = 1:4, 
        method = "cut", 
        itemdescrip = itemlabels$itemnr[1:4])

[1] "Please press Zoom on the Plots window to see the plot"

Code

### also suggested:
# library(RASCHplot) # install first with `devtools::install_github("ERRTG/RASCHplot")`
# CICCplot(PCM(df),
#          which.item = c(1:4),
#          lower.groups = c(0,7,14,21,28,35),
#          grid.items = TRUE)

Code

ICCplot(as.data.frame(df_r), 
        itemnumber = 5:8, 
        method = "cut", 
        itemdescrip = itemlabels$itemnr[5:8])

[1] "Please press Zoom on the Plots window to see the plot"

Code

RIrestscore(df_r)

Item	Observed value	Model expected value	Absolute difference	Adjusted p-value (BH)	Statistical significance level	Location	Relative location
MADRS.S_1	0.76	0.65	0.11	0.000	***	0.27	1.79
MADRS.S_2	0.59	0.66	0.07	0.000	***	-0.88	0.65
MADRS.S_4	0.45	0.64	0.19	0.000	***	1.06	2.59
MADRS.S_5	0.65	0.65	0.00	0.844		-0.20	1.32
MADRS.S_6	0.73	0.65	0.08	0.000	***	-0.30	1.22
MADRS.S_7	0.75	0.64	0.11	0.000	***	0.14	1.67
MADRS.S_8	0.64	0.67	0.03	0.000	***	-0.60	0.93
MADRS.S_9	0.70	0.64	0.06	0.000	***	0.61	2.14

Code

RIbootRestscore(df_r,cpu=n_cores,iterations = 500,samplesize=800) # samplesize=600

Item	Item-restscore result	% of iterations	Conditional MSQ infit	Relative average item location
MADRS.S_1	overfit	100.0	0.72	1.79
MADRS.S_2	no misfit	10.2	1.23	0.65
MADRS.S_2	underfit	89.8	1.23	0.65
MADRS.S_4	underfit	100.0	1.50	2.59
MADRS.S_5	no misfit	98.0	1.03	1.32
MADRS.S_6	overfit	99.4	0.83	1.22
MADRS.S_7	overfit	100.0	0.77	1.67
MADRS.S_8	no misfit	67.8	1.13	0.93
MADRS.S_8	underfit	32.2	1.13	0.93
MADRS.S_9	no misfit	10.6	0.89	2.14
MADRS.S_9	overfit	89.4	0.89	2.14
Note:
Results based on 500 bootstrap iterations with a sample size of 800. Conditional mean-square infit based on complete responders only, n = 5465.

Code

RIpcmPCA(df_r)

PCA of Rasch model residuals
Eigenvalues	Proportion of variance
1.58	23.8%
1.53	18.8%
1.19	17.8%
1.11	12.9%
0.90	9.9%

Code

simcor1 <- RIgetResidCor(df_r, iterations = 500, cpu = n_cores)
RIresidcorr(df_r, cutoff = simcor1$p99)

	MADRS.S_1	MADRS.S_2	MADRS.S_4	MADRS.S_5	MADRS.S_6	MADRS.S_7	MADRS.S_8
MADRS.S_1
MADRS.S_2	-0.1
MADRS.S_4	-0.18	-0.15
MADRS.S_5	-0.21	-0.11	-0.13
MADRS.S_6	-0.13	-0.22	-0.15	0.1
MADRS.S_7	0.09	-0.31	-0.16	-0.11	0.03
MADRS.S_8	-0.16	-0.04	-0.28	-0.26	-0.23	-0.17
MADRS.S_9	0.06	-0.25	-0.13	-0.25	-0.18	0.04	0.02
Note:
Relative cut-off value is -0.066, which is 0.061 above the average correlation (-0.127). Correlations above the cut-off are highlighted in red text.

Code

RIpartgamLD(df_r)

Item 1	Item 2	Partial gamma	SE	Lower CI	Upper CI	Adjusted p-value (BH)
MADRS.S_6	MADRS.S_5	0.343	0.020	0.305	0.382	0.000
MADRS.S_1	MADRS.S_7	0.295	0.022	0.251	0.339	0.000
MADRS.S_1	MADRS.S_9	0.293	0.022	0.249	0.336	0.000
MADRS.S_5	MADRS.S_6	0.283	0.021	0.243	0.323	0.000
MADRS.S_7	MADRS.S_1	0.283	0.022	0.240	0.327	0.000
MADRS.S_7	MADRS.S_9	0.278	0.023	0.233	0.322	0.000
MADRS.S_7	MADRS.S_6	0.249	0.022	0.206	0.293	0.000
MADRS.S_9	MADRS.S_7	0.247	0.023	0.201	0.292	0.000
MADRS.S_9	MADRS.S_1	0.245	0.022	0.201	0.288	0.000
MADRS.S_9	MADRS.S_8	0.239	0.020	0.199	0.280	0.000
MADRS.S_6	MADRS.S_7	0.226	0.023	0.182	0.270	0.000
MADRS.S_8	MADRS.S_9	0.168	0.021	0.127	0.208	0.000
MADRS.S_8	MADRS.S_2	0.162	0.019	0.126	0.199	0.000
MADRS.S_2	MADRS.S_8	0.137	0.019	0.100	0.173	0.000
MADRS.S_1	MADRS.S_2	0.107	0.022	0.064	0.150	0.000
MADRS.S_1	MADRS.S_4	0.096	0.023	0.051	0.141	0.002
MADRS.S_6	MADRS.S_4	0.092	0.022	0.049	0.136	0.002
MADRS.S_7	MADRS.S_4	0.087	0.023	0.041	0.133	0.011
MADRS.S_7	MADRS.S_5	0.083	0.023	0.038	0.127	0.015
MADRS.S_5	MADRS.S_2	0.08	0.020	0.041	0.119	0.003
MADRS.S_1	MADRS.S_8	0.074	0.021	0.033	0.116	0.025

Code

RIloadLoc(df_r)

Code

mirt(df_r, model=1, itemtype='Rasch', verbose = FALSE) %>% 
  plot(type="trace", as.table = TRUE, 
       theta_lim = c(-6,6))

Code

# for fewer items or a more magnified figure, use:
#RIitemCats(df)

Code

# increase fig-height above as needed, if you have many items
RItargeting(df_r)

Code

RIitemHierarchy(df_r)

Code

iarm::score_groups(as.data.frame(df_r)) %>% 
  as.data.frame(nm = "score_group") %>% 
  dplyr::count(score_group)

  score_group    n
1           1 2709
2           2 2756

Code

dif_plots <- df_r %>% 
  add_column(dif = iarm::score_groups(.)) %>% 
  split(.$dif) %>% # split the data using the DIF variable
  map(~ RItileplot(.x %>% dplyr::select(!dif)) + labs(title = .x$dif))
dif_plots[[1]] + dif_plots[[2]]

Code

clr_tests(df_r, model = "PCM")


Conditional Likelihood Ratio Tests:

        clr df pvalue sig 
overall 967 46 0       ***

Code

item_obsexp(PCM(df_r))

Score group 1: 
          mean obs mean exp std.res sig
MADRS.S_1  0.709    0.792   -6.435  -- 
MADRS.S_2  1.826    1.758    4.655  ++ 
MADRS.S_4  0.347    0.209   15.423  ++ 
MADRS.S_5  1.172    1.157    1.107     
MADRS.S_6  1.098    1.145   -3.484  -- 
MADRS.S_7  0.892    0.972   -6.087  -- 
MADRS.S_8  1.401    1.371    1.921     
MADRS.S_9  0.559    0.600   -3.222  -- 

Score group 2: 
          mean obs mean exp std.res sig
MADRS.S_1  2.418    2.340    4.542  ++ 
MADRS.S_2  3.286    3.349   -4.018  -- 
MADRS.S_4  1.244    1.375   -7.717  -- 
MADRS.S_5  2.717    2.731   -0.854     
MADRS.S_6  2.821    2.776    2.573  +  
MADRS.S_7  2.406    2.330    4.886  ++ 
MADRS.S_8  3.373    3.402   -1.700     
MADRS.S_9  1.897    1.858    2.620  +

Code

grouping_based_on_score = score_groups(as.data.frame(df_r))
partgam_DIF(dat.items = as.data.frame(df_r),
            dat.exo = grouping_based_on_score)

       Item                     Var gamma  se pvalue padj.BH sig lower upper
1 MADRS.S_1 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
2 MADRS.S_2 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
3 MADRS.S_4 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
4 MADRS.S_5 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
5 MADRS.S_6 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
6 MADRS.S_7 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
7 MADRS.S_8 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
8 MADRS.S_9 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA

Code

dif.sex <- dif %>% filter(Time=='WEEK02') %>% select(sex) %>% as.vector(.)
dif.sex <- dif.sex[[1]] #female = 1, male = 2 
RIdifTable(df_r, dif.sex)

Item	2	3	Mean location	StDev	MaxDiff
MADRS.S_1	0.270	0.226	0.248	0.031	0.044
MADRS.S_2	-0.934	-0.832	-0.883	0.072	0.102
MADRS.S_4	0.974	1.240	1.107	0.188	0.266
MADRS.S_5	-0.183	-0.308	-0.246	0.089	0.125
MADRS.S_6	-0.225	-0.493	-0.359	0.189	0.268
MADRS.S_7	0.196	-0.015	0.090	0.149	0.211
MADRS.S_8	-0.668	-0.521	-0.594	0.104	0.147
MADRS.S_9	0.571	0.704	0.637	0.094	0.133

Code

partgam_DIF(dat.items = as.data.frame(df_r),
            dat.exo = dif.sex)

       Item     Var   gamma     se pvalue padj.BH  sig   lower   upper
1 MADRS.S_1 dif.sex  0.2209 0.0265 0.0000  0.0000  ***  0.1689  0.2728
2 MADRS.S_2 dif.sex -0.1625 0.0239 0.0000  0.0000  *** -0.2093 -0.1156
3 MADRS.S_4 dif.sex -0.1092 0.0273 0.0001  0.0005  *** -0.1626 -0.0557
4 MADRS.S_5 dif.sex  0.0543 0.0255 0.0330  0.2641       0.0044  0.1042
5 MADRS.S_6 dif.sex  0.0838 0.0263 0.0014  0.0116    *  0.0322  0.1354
6 MADRS.S_7 dif.sex  0.2262 0.0269 0.0000  0.0000  ***  0.1735  0.2788
7 MADRS.S_8 dif.sex -0.1769 0.0236 0.0000  0.0000  *** -0.2232 -0.1306
8 MADRS.S_9 dif.sex  0.1004 0.0280 0.0003  0.0026   **  0.0456  0.1552

Code

dif.age <- dif %>% filter(Time=='WEEK02') %>% select(age) %>% as.vector(.)
dif.age <- dif.age[[1]]
RIdifTable(df_r, dif.age)

Item	3	4	6	8	9	Mean location	StDev	MaxDiff
MADRS.S_1	0.281	0.114	-0.216	0.265	-0.001	0.088	0.206	0.497
MADRS.S_2	-0.966	-1.008	-0.862	-0.672	-0.773	-0.856	0.138	0.336
MADRS.S_4	0.883	1.040	1.075	1.004	0.745	0.950	0.135	0.330
MADRS.S_5	-0.228	-0.028	-0.205	-0.248	0.228	-0.096	0.201	0.476
MADRS.S_6	-0.332	-0.330	-0.203	-0.332	-0.373	-0.314	0.065	0.170
MADRS.S_7	0.123	0.382	0.193	-0.010	-0.396	0.058	0.291	0.778
MADRS.S_8	-0.510	-0.740	-0.474	-0.617	-0.161	-0.500	0.216	0.578
MADRS.S_9	0.749	0.570	0.692	0.611	0.731	0.671	0.078	0.180

Code

partgam_DIF(dat.items = as.data.frame(df_r),
            dat.exo = dif.age)

       Item     Var   gamma     se pvalue padj.BH  sig   lower   upper
1 MADRS.S_1 dif.age -0.0067 0.0162 0.6794  1.0000      -0.0384  0.0250
2 MADRS.S_2 dif.age -0.1251 0.0140 0.0000  0.0000  *** -0.1526 -0.0976
3 MADRS.S_4 dif.age -0.0255 0.0156 0.1038  0.8303      -0.0561  0.0052
4 MADRS.S_5 dif.age  0.0386 0.0149 0.0095  0.0760    .  0.0094  0.0677
5 MADRS.S_6 dif.age -0.0020 0.0154 0.8975  1.0000      -0.0321  0.0281
6 MADRS.S_7 dif.age  0.1424 0.0159 0.0000  0.0000  ***  0.1113  0.1736
7 MADRS.S_8 dif.age -0.0070 0.0139 0.6141  1.0000      -0.0342  0.0202
8 MADRS.S_9 dif.age  0.0941 0.0163 0.0000  0.0000  ***  0.0622  0.1260

Code

dif.edu <- dif %>% filter(Time=='WEEK02') %>% select(education) %>% as.vector(.)
dif.edu <- dif.edu[[1]]
RIdifTable(df_r, dif.edu)

Item	2	3	Mean location	StDev	MaxDiff
MADRS.S_1	0.272	0.213	0.242	0.041	0.058
MADRS.S_2	-0.890	-0.887	-0.889	0.002	0.003
MADRS.S_4	0.951	1.241	1.096	0.205	0.290
MADRS.S_5	-0.241	-0.187	-0.214	0.038	0.054
MADRS.S_6	-0.374	-0.222	-0.298	0.107	0.152
MADRS.S_7	0.144	0.102	0.123	0.030	0.042
MADRS.S_8	-0.522	-0.761	-0.641	0.169	0.239
MADRS.S_9	0.661	0.502	0.581	0.112	0.159

Code

partgam_DIF(dat.items = as.data.frame(df_r),
            dat.exo = dif.edu)

       Item     Var   gamma     se pvalue padj.BH  sig   lower   upper
1 MADRS.S_1 dif.edu  0.0287 0.0243 0.2381  1.0000      -0.0190  0.0763
2 MADRS.S_2 dif.edu -0.0383 0.0214 0.0738  0.5902      -0.0802  0.0037
3 MADRS.S_4 dif.edu -0.0618 0.0238 0.0095  0.0761    . -0.1085 -0.0151
4 MADRS.S_5 dif.edu -0.0232 0.0229 0.3096  1.0000      -0.0681  0.0216
5 MADRS.S_6 dif.edu -0.0912 0.0235 0.0001  0.0008  *** -0.1371 -0.0452
6 MADRS.S_7 dif.edu  0.0730 0.0247 0.0031  0.0246    *  0.0247  0.1214
7 MADRS.S_8 dif.edu  0.0894 0.0216 0.0000  0.0003  ***  0.0471  0.1317
8 MADRS.S_9 dif.edu  0.0582 0.0250 0.0198  0.1580       0.0093  0.1071

Code

dif.yr <- dif %>% filter(Time=='WEEK02') %>% select(TreatmentAccessStart) %>% as.vector(.)
dif.yr <- dif.yr[[1]]
RIdifTable(df_r, dif.yr)

[1] "No statistically significant DIF found."

Code

partgam_DIF(dat.items = as.data.frame(df_r),
            dat.exo = dif.yr)

       Item    Var   gamma     se pvalue padj.BH  sig   lower   upper
1 MADRS.S_1 dif.yr -0.0159 0.0167 0.3418  1.0000      -0.0486  0.0168
2 MADRS.S_2 dif.yr  0.0562 0.0150 0.0002  0.0015   **  0.0267  0.0856
3 MADRS.S_4 dif.yr -0.0241 0.0168 0.1511  1.0000      -0.0569  0.0088
4 MADRS.S_5 dif.yr  0.0086 0.0157 0.5856  1.0000      -0.0222  0.0394
5 MADRS.S_6 dif.yr -0.0024 0.0163 0.8849  1.0000      -0.0342  0.0295
6 MADRS.S_7 dif.yr -0.0183 0.0170 0.2818  1.0000      -0.0517  0.0151
7 MADRS.S_8 dif.yr  0.0219 0.0149 0.1399  1.0000      -0.0072  0.0510
8 MADRS.S_9 dif.yr -0.0754 0.0170 0.0000  0.0001  *** -0.1087 -0.0421

9.1 Analysis part 2 decision

Remove item 4, appetite : from the residual contrast loadings we see it loads to something completely different than the others.

10 Rasch analysis 3

Removing item 4

Code

df_r3 <- df_week02 %>% select(!c(MADRS.S_4,MADRS.S_3))

Code

RIlistItemsMargin(df_r3, fontsize = 12)

itemnr	item
MADRS.S_1	Sadness
MADRS.S_2	Anxiety
MADRS.S_5	Concentration difficulties
MADRS.S_6	Lassitude
MADRS.S_7	Inability to feel
MADRS.S_8	Pessimism
MADRS.S_9	Suicidal thoughts

Code

RIitemfit(df_r3, cutoff = "Smith98")

Item	InfitMSQ	Location	1 ± 2 / √n
MADRS.S_1	0.775	0.45	[0.973, 1.027]
MADRS.S_2	1.325	-0.77	[0.973, 1.027]
MADRS.S_5	1.119	-0.05	[0.973, 1.027]
MADRS.S_6	0.901	-0.15	[0.973, 1.027]
MADRS.S_7	0.823	0.34	[0.973, 1.027]
MADRS.S_8	1.169	-0.49	[0.973, 1.027]
MADRS.S_9	0.958	0.81	[0.973, 1.027]
Note:
MSQ values based on conditional estimation (n = 5465 complete cases).

Code

simfit1 <- RIgetfit(df_r3, iterations = 300, cpu = n_cores)  

RIitemfit(df_r3, simfit1)

Item	InfitMSQ	Infit thresholds	OutfitMSQ	Outfit thresholds	Infit diff	Outfit diff	Location
MADRS.S_1	0.775	[0.952, 1.049]	0.76	[0.949, 1.06]	0.177	0.189	0.45
MADRS.S_2	1.325	[0.963, 1.038]	1.342	[0.959, 1.041]	0.287	0.301	-0.77
MADRS.S_5	1.119	[0.955, 1.047]	1.105	[0.956, 1.046]	0.072	0.059	-0.05
MADRS.S_6	0.901	[0.952, 1.052]	0.896	[0.954, 1.046]	0.051	0.058	-0.15
MADRS.S_7	0.823	[0.954, 1.043]	0.807	[0.952, 1.046]	0.131	0.145	0.34
MADRS.S_8	1.169	[0.95, 1.046]	1.194	[0.95, 1.048]	0.123	0.146	-0.49
MADRS.S_9	0.958	[0.955, 1.051]	0.908	[0.938, 1.073]	no misfit	0.030	0.81
Note:
MSQ values based on conditional calculations (n = 5465 complete cases). Simulation based thresholds from 300 simulated datasets.

Code

RIgetfitPlot(simfit1, df_r3)

Code

ICCplot(as.data.frame(df_r3), 
        itemnumber = 1:4, 
        method = "cut", 
        itemdescrip = itemlabels$itemnr[1:4])

[1] "Please press Zoom on the Plots window to see the plot"

Code

### also suggested:
# library(RASCHplot) # install first with `devtools::install_github("ERRTG/RASCHplot")`
# CICCplot(PCM(df),
#          which.item = c(1:4),
#          lower.groups = c(0,7,14,21,28,35),
#          grid.items = TRUE)

Code

ICCplot(as.data.frame(df_r3), 
        itemnumber = 5:7, 
        method = "cut", 
        itemdescrip = itemlabels$itemnr[5:7])

[1] "Please press Zoom on the Plots window to see the plot"

Code

RIrestscore(df_r3)

Item	Observed value	Model expected value	Absolute difference	Adjusted p-value (BH)	Statistical significance level	Location	Relative location
MADRS.S_1	0.76	0.67	0.09	0.000	***	0.45	1.94
MADRS.S_2	0.59	0.68	0.09	0.000	***	-0.77	0.71
MADRS.S_5	0.65	0.68	0.03	0.002	**	-0.05	1.44
MADRS.S_6	0.72	0.68	0.04	0.000	***	-0.15	1.33
MADRS.S_7	0.75	0.67	0.08	0.000	***	0.34	1.82
MADRS.S_8	0.65	0.70	0.05	0.000	***	-0.49	0.99
MADRS.S_9	0.70	0.66	0.04	0.000	***	0.81	2.29

Code

RIbootRestscore(df_r3,cpu=n_cores,iterations = 500,samplesize=800) # samplesize=600

Item	Item-restscore result	% of iterations	Conditional MSQ infit	Relative average item location
MADRS.S_1	overfit	99.8	0.78	1.94
MADRS.S_2	underfit	99.8	1.33	0.71
MADRS.S_5	no misfit	90.0	1.12	1.44
MADRS.S_5	underfit	10.0	1.12	1.44
MADRS.S_6	no misfit	25.4	0.90	1.33
MADRS.S_6	overfit	74.6	0.90	1.33
MADRS.S_7	overfit	99.8	0.82	1.82
MADRS.S_8	no misfit	52.0	1.17	0.99
MADRS.S_8	underfit	48.0	1.17	0.99
MADRS.S_9	no misfit	62.0	0.96	2.29
MADRS.S_9	overfit	38.0	0.96	2.29
Note:
Results based on 500 bootstrap iterations with a sample size of 800. Conditional mean-square infit based on complete responders only, n = 5465.

Code

RIpcmPCA(df_r3)

PCA of Rasch model residuals
Eigenvalues	Proportion of variance
1.61	24.7%
1.57	23.5%
1.18	17%
0.93	12.9%
0.88	11.4%

Code

simcor1 <- RIgetResidCor(df_r3, iterations = 500, cpu = n_cores)
RIresidcorr(df_r3, cutoff = simcor1$p99)

	MADRS.S_1	MADRS.S_2	MADRS.S_5	MADRS.S_6	MADRS.S_7	MADRS.S_8
MADRS.S_1
MADRS.S_2	-0.13
MADRS.S_5	-0.25	-0.13
MADRS.S_6	-0.17	-0.25	0.08
MADRS.S_7	0.06	-0.34	-0.14	0
MADRS.S_8	-0.22	-0.08	-0.31	-0.28	-0.22
MADRS.S_9	0.03	-0.28	-0.28	-0.21	0.01	-0.02
Note:
Relative cut-off value is -0.086, which is 0.063 above the average correlation (-0.149). Correlations above the cut-off are highlighted in red text.

Code

RIpartgamLD(df_r3)

Item 1	Item 2	Partial gamma	SE	Lower CI	Upper CI	Adjusted p-value (BH)
MADRS.S_6	MADRS.S_5	0.346	0.020	0.308	0.385	0.000
MADRS.S_1	MADRS.S_7	0.298	0.022	0.254	0.341	0.000
MADRS.S_1	MADRS.S_9	0.295	0.022	0.252	0.338	0.000
MADRS.S_7	MADRS.S_1	0.284	0.022	0.240	0.327	0.000
MADRS.S_5	MADRS.S_6	0.281	0.020	0.241	0.321	0.000
MADRS.S_7	MADRS.S_9	0.275	0.023	0.230	0.319	0.000
MADRS.S_7	MADRS.S_6	0.25	0.022	0.207	0.293	0.000
MADRS.S_9	MADRS.S_1	0.243	0.022	0.199	0.286	0.000
MADRS.S_9	MADRS.S_7	0.242	0.023	0.197	0.287	0.000
MADRS.S_6	MADRS.S_7	0.23	0.022	0.186	0.274	0.000
MADRS.S_9	MADRS.S_8	0.208	0.021	0.167	0.248	0.000
MADRS.S_8	MADRS.S_9	0.141	0.021	0.100	0.183	0.000
MADRS.S_8	MADRS.S_2	0.138	0.019	0.101	0.176	0.000
MADRS.S_2	MADRS.S_8	0.101	0.019	0.064	0.139	0.000
MADRS.S_1	MADRS.S_2	0.099	0.022	0.056	0.142	0.000
MADRS.S_7	MADRS.S_5	0.084	0.023	0.039	0.128	0.009
MADRS.S_5	MADRS.S_2	0.07	0.020	0.031	0.109	0.021

Code

RIloadLoc(df_r3)

Code

mirt(df_r3, model=1, itemtype='Rasch', verbose = FALSE) %>% 
  plot(type="trace", as.table = TRUE, 
       theta_lim = c(-6,6))

Code

# for fewer items or a more magnified figure, use:
#RIitemCats(df)

Code

# increase fig-height above as needed, if you have many items
RItargeting(df_r3)

Code

RIitemHierarchy(df_r3)

Code

iarm::score_groups(as.data.frame(df_r3)) %>% 
  as.data.frame(nm = "score_group") %>% 
  dplyr::count(score_group)

  score_group    n
1           1 2898
2           2 2567

Code

dif_plots <- df_r3 %>% 
  add_column(dif = iarm::score_groups(.)) %>% 
  split(.$dif) %>% # split the data using the DIF variable
  map(~ RItileplot(.x %>% dplyr::select(!dif)) + labs(title = .x$dif))
dif_plots[[1]] + dif_plots[[2]]

Code

clr_tests(df_r3, model = "PCM")


Conditional Likelihood Ratio Tests:

        clr df pvalue sig 
overall 595 40 0       ***

Code

item_obsexp(PCM(df_r3))

Score group 1: 
          mean obs mean exp std.res sig
MADRS.S_1  0.755    0.817   -5.036  -- 
MADRS.S_2  1.866    1.792    5.391  ++ 
MADRS.S_5  1.223    1.187    2.773  +  
MADRS.S_6  1.152    1.173   -1.635     
MADRS.S_7  0.945    0.999   -4.325  -- 
MADRS.S_8  1.459    1.418    2.722  +  
MADRS.S_9  0.608    0.621   -1.092     

Score group 2: 
          mean obs mean exp std.res sig
MADRS.S_1  2.495    2.428    3.940  ++ 
MADRS.S_2  3.353    3.433   -5.121  -- 
MADRS.S_5  2.777    2.816   -2.356  -  
MADRS.S_6  2.891    2.868    1.343     
MADRS.S_7  2.462    2.403    3.823  ++ 
MADRS.S_8  3.458    3.502   -2.705  -  
MADRS.S_9  1.943    1.929    0.989

Code

grouping_based_on_score = score_groups(as.data.frame(df_r3))
partgam_DIF(dat.items = as.data.frame(df_r3),
            dat.exo = grouping_based_on_score)

       Item                     Var gamma  se pvalue padj.BH sig lower upper
1 MADRS.S_1 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
2 MADRS.S_2 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
3 MADRS.S_5 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
4 MADRS.S_6 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
5 MADRS.S_7 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
6 MADRS.S_8 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
7 MADRS.S_9 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA

Code

dif.sex <- dif %>% filter(Time=='WEEK02') %>% select(sex) %>% as.vector(.)
dif.sex <- dif.sex[[1]] #female = 1, male = 2 
RIdifTable(df_r3, dif.sex)

Item	2	3	Mean location	StDev	MaxDiff
MADRS.S_1	0.443	0.423	0.433	0.014	0.020
MADRS.S_2	-0.855	-0.693	-0.774	0.114	0.162
MADRS.S_5	-0.044	-0.129	-0.087	0.060	0.085
MADRS.S_6	-0.089	-0.325	-0.207	0.166	0.235
MADRS.S_7	0.375	0.191	0.283	0.130	0.184
MADRS.S_8	-0.584	-0.377	-0.480	0.146	0.207
MADRS.S_9	0.754	0.909	0.831	0.110	0.155

Code

partgam_DIF(dat.items = as.data.frame(df_r3),
            dat.exo = dif.sex)

       Item     Var   gamma     se pvalue padj.BH  sig   lower   upper
1 MADRS.S_1 dif.sex  0.1993 0.0267 0.0000  0.0000  ***  0.1469  0.2517
2 MADRS.S_2 dif.sex -0.1754 0.0239 0.0000  0.0000  *** -0.2222 -0.1286
3 MADRS.S_5 dif.sex  0.0365 0.0255 0.1526  1.0000      -0.0135  0.0865
4 MADRS.S_6 dif.sex  0.0685 0.0264 0.0095  0.0667    .  0.0167  0.1202
5 MADRS.S_7 dif.sex  0.2089 0.0271 0.0000  0.0000  ***  0.1558  0.2619
6 MADRS.S_8 dif.sex -0.2011 0.0237 0.0000  0.0000  *** -0.2476 -0.1547
7 MADRS.S_9 dif.sex  0.0903 0.0280 0.0013  0.0088   **  0.0354  0.1452

Code

dif.age <- dif %>% filter(Time=='WEEK02') %>% select(age) %>% as.vector(.)
dif.age <- dif.age[[1]]
RIdifTable(df_r3, dif.age)

Item	3	4	6	8	9	Mean location	StDev	MaxDiff
MADRS.S_1	0.428	0.301	-0.098	0.449	0.016	0.219	0.248	0.547
MADRS.S_2	-0.885	-0.915	-0.784	-0.479	-0.885	-0.789	0.181	0.436
MADRS.S_5	-0.101	0.120	-0.018	-0.139	0.255	0.023	0.163	0.394
MADRS.S_6	-0.208	-0.183	-0.026	-0.233	-0.428	-0.216	0.144	0.403
MADRS.S_7	0.277	0.563	0.408	0.146	0.384	0.356	0.156	0.417
MADRS.S_8	-0.414	-0.648	-0.405	-0.537	-0.102	-0.421	0.204	0.546
MADRS.S_9	0.905	0.760	0.923	0.794	0.761	0.828	0.079	0.163

Code

partgam_DIF(dat.items = as.data.frame(df_r3),
            dat.exo = dif.age)

       Item     Var   gamma     se pvalue padj.BH  sig   lower   upper
1 MADRS.S_1 dif.age -0.0134 0.0162 0.4054  1.0000      -0.0451  0.0182
2 MADRS.S_2 dif.age -0.1291 0.0140 0.0000  0.0000  *** -0.1565 -0.1016
3 MADRS.S_5 dif.age  0.0339 0.0148 0.0222  0.1555       0.0049  0.0630
4 MADRS.S_6 dif.age -0.0081 0.0154 0.5984  1.0000      -0.0384  0.0221
5 MADRS.S_7 dif.age  0.1370 0.0159 0.0000  0.0000  ***  0.1058  0.1682
6 MADRS.S_8 dif.age -0.0097 0.0140 0.4882  1.0000      -0.0371  0.0177
7 MADRS.S_9 dif.age  0.0904 0.0163 0.0000  0.0000  ***  0.0584  0.1223

Code

dif.edu <- dif %>% filter(Time=='WEEK02') %>% select(education) %>% as.vector(.)
dif.edu <- dif.edu[[1]]
RIdifTable(df_r3, dif.edu)

Item	2	3	Mean location	StDev	MaxDiff
MADRS.S_1	0.433	0.427	0.430	0.004	0.006
MADRS.S_2	-0.802	-0.768	-0.785	0.024	0.034
MADRS.S_5	-0.107	-0.008	-0.057	0.070	0.099
MADRS.S_6	-0.248	-0.043	-0.146	0.145	0.205
MADRS.S_7	0.312	0.320	0.316	0.006	0.008
MADRS.S_8	-0.424	-0.648	-0.536	0.159	0.225
MADRS.S_9	0.835	0.720	0.778	0.082	0.115

Code

partgam_DIF(dat.items = as.data.frame(df_r3),
            dat.exo = dif.edu)

       Item     Var   gamma     se pvalue padj.BH  sig   lower   upper
1 MADRS.S_1 dif.edu  0.0136 0.0243 0.5752  1.0000      -0.0340  0.0613
2 MADRS.S_2 dif.edu -0.0427 0.0214 0.0464  0.3247      -0.0848 -0.0007
3 MADRS.S_5 dif.edu -0.0330 0.0228 0.1478  1.0000      -0.0777  0.0117
4 MADRS.S_6 dif.edu -0.1015 0.0234 0.0000  0.0001  *** -0.1474 -0.0556
5 MADRS.S_7 dif.edu  0.0671 0.0246 0.0063  0.0439    *  0.0190  0.1153
6 MADRS.S_8 dif.edu  0.0909 0.0218 0.0000  0.0002  ***  0.0482  0.1336
7 MADRS.S_9 dif.edu  0.0493 0.0249 0.0477  0.3338       0.0005  0.0982

Code

dif.yr <- dif %>% filter(Time=='WEEK02') %>% select(TreatmentAccessStart) %>% as.vector(.)
dif.yr <- dif.yr[[1]]
RIdifTable(df_r3, dif.yr)

[1] "No statistically significant DIF found."

Code

partgam_DIF(dat.items = as.data.frame(df_r3),
            dat.exo = dif.yr)

       Item    Var   gamma     se pvalue padj.BH  sig   lower   upper
1 MADRS.S_1 dif.yr -0.0195 0.0167 0.2413  1.0000      -0.0522  0.0131
2 MADRS.S_2 dif.yr  0.0509 0.0151 0.0007  0.0052   **  0.0213  0.0805
3 MADRS.S_5 dif.yr  0.0069 0.0157 0.6613  1.0000      -0.0239  0.0376
4 MADRS.S_6 dif.yr -0.0063 0.0163 0.7005  1.0000      -0.0383  0.0257
5 MADRS.S_7 dif.yr -0.0182 0.0170 0.2862  1.0000      -0.0515  0.0152
6 MADRS.S_8 dif.yr  0.0190 0.0151 0.2063  1.0000      -0.0105  0.0486
7 MADRS.S_9 dif.yr -0.0780 0.0170 0.0000  0.0000  *** -0.1113 -0.0448

10.1 Analysis part 3 decision

High number of boot misfit for item 2 and for conditional item fit it was most deviant - removing.

11 Rasch analysis 4

Removing item 2

Code

df_r3 <- df_week02 %>% select(!c(MADRS.S_4,MADRS.S_3,MADRS.S_2))

Code

RIlistItemsMargin(df_r3, fontsize = 12)

itemnr	item
MADRS.S_1	Sadness
MADRS.S_5	Concentration difficulties
MADRS.S_6	Lassitude
MADRS.S_7	Inability to feel
MADRS.S_8	Pessimism
MADRS.S_9	Suicidal thoughts

Code

RIitemfit(df_r3, cutoff = "Smith98")

Item	InfitMSQ	Location	1 ± 2 / √n
MADRS.S_1	0.838	0.34	[0.973, 1.027]
MADRS.S_5	1.215	-0.19	[0.973, 1.027]
MADRS.S_6	0.929	-0.30	[0.973, 1.027]
MADRS.S_7	0.806	0.22	[0.973, 1.027]
MADRS.S_8	1.321	-0.66	[0.973, 1.027]
MADRS.S_9	0.967	0.70	[0.973, 1.027]
Note:
MSQ values based on conditional estimation (n = 5465 complete cases).

Code

simfit1 <- RIgetfit(df_r3, iterations = 300, cpu = n_cores)  

RIitemfit(df_r3, simfit1)

Item	InfitMSQ	Infit thresholds	OutfitMSQ	Outfit thresholds	Infit diff	Outfit diff	Location
MADRS.S_1	0.838	[0.95, 1.05]	0.843	[0.943, 1.054]	0.112	0.100	0.34
MADRS.S_5	1.215	[0.949, 1.059]	1.214	[0.947, 1.057]	0.156	0.157	-0.19
MADRS.S_6	0.929	[0.953, 1.04]	0.93	[0.954, 1.038]	0.024	0.024	-0.30
MADRS.S_7	0.806	[0.955, 1.048]	0.793	[0.955, 1.05]	0.149	0.162	0.22
MADRS.S_8	1.321	[0.955, 1.057]	1.358	[0.954, 1.057]	0.264	0.301	-0.66
MADRS.S_9	0.967	[0.959, 1.049]	0.923	[0.951, 1.067]	no misfit	0.028	0.70
Note:
MSQ values based on conditional calculations (n = 5465 complete cases). Simulation based thresholds from 300 simulated datasets.

Code

RIgetfitPlot(simfit1, df_r3)

Code

ICCplot(as.data.frame(df_r3), 
        itemnumber = 1:4, 
        method = "cut", 
        itemdescrip = itemlabels$itemnr[1:4])

[1] "Please press Zoom on the Plots window to see the plot"

Code

### also suggested:
# library(RASCHplot) # install first with `devtools::install_github("ERRTG/RASCHplot")`
# CICCplot(PCM(df),
#          which.item = c(1:4),
#          lower.groups = c(0,7,14,21,28,35),
#          grid.items = TRUE)

Code

ICCplot(as.data.frame(df_r3), 
        itemnumber = 5:6, 
        method = "cut", 
        itemdescrip = itemlabels$itemnr[5:6])

[1] "Please press Zoom on the Plots window to see the plot"

Code

RIrestscore(df_r3)

Item	Observed value	Model expected value	Absolute difference	Statistical significance level	Location	Relative location
MADRS.S_1	0.76	0.69	0.07	***	0.34	2.03
MADRS.S_5	0.64	0.70	0.06	***	-0.19	1.50
MADRS.S_6	0.73	0.70	0.03	***	-0.30	1.40
MADRS.S_7	0.77	0.69	0.08	***	0.22	1.91
MADRS.S_8	0.64	0.71	0.07	***	-0.66	1.03
MADRS.S_9	0.72	0.69	0.03	***	0.70	2.40

Code

RIbootRestscore(df_r3,cpu=n_cores,iterations = 500,samplesize=800) # samplesize=600

Item	Item-restscore result	% of iterations	Conditional MSQ infit	Relative average item location
MADRS.S_1	overfit	96.0	0.84	2.03
MADRS.S_5	no misfit	25.6	1.21	1.50
MADRS.S_5	underfit	74.4	1.21	1.50
MADRS.S_6	no misfit	53.4	0.93	1.40
MADRS.S_6	overfit	46.6	0.93	1.40
MADRS.S_7	overfit	100.0	0.81	1.91
MADRS.S_8	underfit	98.6	1.32	1.03
MADRS.S_9	no misfit	68.8	0.97	2.40
MADRS.S_9	overfit	31.2	0.97	2.40
Note:
Results based on 500 bootstrap iterations with a sample size of 800. Conditional mean-square infit based on complete responders only, n = 5465.

Code

RIpcmPCA(df_r3)

PCA of Rasch model residuals
Eigenvalues	Proportion of variance
1.66	30.8%
1.37	23.5%
1.03	16.7%
0.97	15.1%
0.95	13.7%

Code

simcor1 <- RIgetResidCor(df_r3, iterations = 500, cpu = n_cores)
RIresidcorr(df_r3, cutoff = simcor1$p99)

	MADRS.S_1	MADRS.S_5	MADRS.S_6	MADRS.S_7	MADRS.S_8
MADRS.S_1
MADRS.S_5	-0.27
MADRS.S_6	-0.22	0.04
MADRS.S_7	0	-0.2	-0.09
MADRS.S_8	-0.23	-0.32	-0.32	-0.28
MADRS.S_9	-0.02	-0.34	-0.29	-0.07	-0.06
Note:
Relative cut-off value is -0.114, which is 0.064 above the average correlation (-0.178). Correlations above the cut-off are highlighted in red text.

Code

RIpartgamLD(df_r3)

Item 1	Item 2	Partial gamma	SE	Lower CI	Upper CI
MADRS.S_6	MADRS.S_5	0.363	0.020	0.324	0.401
MADRS.S_5	MADRS.S_6	0.266	0.021	0.226	0.306
MADRS.S_7	MADRS.S_1	0.258	0.023	0.213	0.303
MADRS.S_1	MADRS.S_9	0.247	0.023	0.202	0.292
MADRS.S_1	MADRS.S_7	0.239	0.023	0.194	0.285
MADRS.S_9	MADRS.S_8	0.222	0.021	0.180	0.263
MADRS.S_9	MADRS.S_1	0.212	0.023	0.167	0.257
MADRS.S_7	MADRS.S_6	0.19	0.023	0.144	0.236
MADRS.S_7	MADRS.S_9	0.186	0.025	0.138	0.235
MADRS.S_9	MADRS.S_7	0.135	0.025	0.087	0.184
MADRS.S_6	MADRS.S_7	0.133	0.024	0.087	0.180
MADRS.S_8	MADRS.S_9	0.1	0.022	0.057	0.142

Code

RIloadLoc(df_r3)

Code

mirt(df_r3, model=1, itemtype='Rasch', verbose = FALSE) %>% 
  plot(type="trace", as.table = TRUE, 
       theta_lim = c(-6,6))

Code

# for fewer items or a more magnified figure, use:
#RIitemCats(df)

Code

# increase fig-height above as needed, if you have many items
RItargeting(df_r3)

Code

RIitemHierarchy(df_r3)

Code

iarm::score_groups(as.data.frame(df_r3)) %>% 
  as.data.frame(nm = "score_group") %>% 
  dplyr::count(score_group)

  score_group    n
1           1 2798
2           2 2667

Code

dif_plots <- df_r3 %>% 
  add_column(dif = iarm::score_groups(.)) %>% 
  split(.$dif) %>% # split the data using the DIF variable
  map(~ RItileplot(.x %>% dplyr::select(!dif)) + labs(title = .x$dif))
dif_plots[[1]] + dif_plots[[2]]

Code

clr_tests(df_r3, model = "PCM")


Conditional Likelihood Ratio Tests:

        clr df pvalue sig 
overall 523 34 0       ***

Code

item_obsexp(PCM(df_r3))

Score group 1: 
          mean obs mean exp std.res sig
MADRS.S_1  0.744    0.799   -4.514  -- 
MADRS.S_5  1.234    1.181    4.225  ++ 
MADRS.S_6  1.150    1.166   -1.252     
MADRS.S_7  0.937    0.988   -4.135  -- 
MADRS.S_8  1.478    1.399    5.501  ++ 
MADRS.S_9  0.583    0.595   -0.968     

Score group 2: 
          mean obs mean exp std.res sig
MADRS.S_1  2.466    2.413    3.336  ++ 
MADRS.S_5  2.750    2.802   -3.360  -- 
MADRS.S_6  2.867    2.852    0.952     
MADRS.S_7  2.445    2.395    3.441  ++ 
MADRS.S_8  3.414    3.492   -5.020  -- 
MADRS.S_9  1.938    1.927    0.826

Code

grouping_based_on_score = score_groups(as.data.frame(df_r3))
partgam_DIF(dat.items = as.data.frame(df_r3),
            dat.exo = grouping_based_on_score)

       Item                     Var gamma  se pvalue padj.BH sig lower upper
1 MADRS.S_1 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
2 MADRS.S_5 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
3 MADRS.S_6 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
4 MADRS.S_7 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
5 MADRS.S_8 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
6 MADRS.S_9 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA

Code

dif.sex <- dif %>% filter(Time=='WEEK02') %>% select(sex) %>% as.vector(.)
dif.sex <- dif.sex[[1]] #female = 1, male = 2 
RIdifTable(df_r3, dif.sex)

Item	2	3	Mean location	StDev	MaxDiff
MADRS.S_1	0.323	0.324	0.324	0.001	0.001
MADRS.S_5	-0.196	-0.255	-0.225	0.042	0.059
MADRS.S_6	-0.243	-0.457	-0.350	0.152	0.214
MADRS.S_7	0.246	0.086	0.166	0.114	0.161
MADRS.S_8	-0.770	-0.520	-0.645	0.177	0.251
MADRS.S_9	0.639	0.821	0.730	0.129	0.182

Code

partgam_DIF(dat.items = as.data.frame(df_r3),
            dat.exo = dif.sex)

       Item     Var   gamma     se pvalue padj.BH  sig   lower   upper
1 MADRS.S_1 dif.sex  0.1585 0.0272 0.0000  0.0000  ***  0.1052  0.2117
2 MADRS.S_5 dif.sex  0.0048 0.0256 0.8518  1.0000      -0.0454  0.0549
3 MADRS.S_6 dif.sex  0.0207 0.0270 0.4420  1.0000      -0.0321  0.0736
4 MADRS.S_7 dif.sex  0.1807 0.0282 0.0000  0.0000  ***  0.1254  0.2360
5 MADRS.S_8 dif.sex -0.2368 0.0234 0.0000  0.0000  *** -0.2826 -0.1909
6 MADRS.S_9 dif.sex  0.0619 0.0287 0.0310  0.1859       0.0057  0.1181

Code

dif.age <- dif %>% filter(Time=='WEEK02') %>% select(age) %>% as.vector(.)
dif.age <- dif.age[[1]]
RIdifTable(df_r3, dif.age)

Item	3	4	6	7	Mean location	StDev	MaxDiff
MADRS.S_1	0.289	0.161	0.456	-0.188	0.180	0.273	0.644
MADRS.S_5	-0.251	-0.033	-0.273	-0.006	-0.141	0.141	0.267
MADRS.S_6	-0.362	-0.348	-0.306	-0.045	-0.265	0.149	0.317
MADRS.S_7	0.137	0.427	0.100	0.104	0.192	0.157	0.327
MADRS.S_8	-0.581	-0.839	-0.703	-0.442	-0.641	0.170	0.397
MADRS.S_9	0.769	0.632	0.726	0.577	0.676	0.087	0.192

Code

partgam_DIF(dat.items = as.data.frame(df_r3),
            dat.exo = dif.age)

       Item     Var   gamma     se pvalue padj.BH  sig   lower   upper
1 MADRS.S_1 dif.age -0.0446 0.0162 0.0059  0.0355    * -0.0764 -0.0128
2 MADRS.S_5 dif.age  0.0128 0.0148 0.3903  1.0000      -0.0163  0.0419
3 MADRS.S_6 dif.age -0.0379 0.0156 0.0153  0.0920    . -0.0685 -0.0073
4 MADRS.S_7 dif.age  0.1152 0.0165 0.0000  0.0000  ***  0.0828  0.1475
5 MADRS.S_8 dif.age -0.0370 0.0140 0.0082  0.0493    * -0.0645 -0.0096
6 MADRS.S_9 dif.age  0.0624 0.0167 0.0002  0.0012   **  0.0296  0.0952

Code

dif.edu <- dif %>% filter(Time=='WEEK02') %>% select(education) %>% as.vector(.)
dif.edu <- dif.edu[[1]]
RIdifTable(df_r3, dif.edu)

Item	2	3	Mean location	StDev	MaxDiff
MADRS.S_1	0.318	0.323	0.321	0.003	0.005
MADRS.S_5	-0.251	-0.144	-0.198	0.076	0.108
MADRS.S_6	-0.397	-0.181	-0.289	0.153	0.216
MADRS.S_7	0.191	0.207	0.199	0.011	0.016
MADRS.S_8	-0.590	-0.827	-0.708	0.167	0.237
MADRS.S_9	0.728	0.621	0.675	0.076	0.108

Code

partgam_DIF(dat.items = as.data.frame(df_r3),
            dat.exo = dif.edu)

       Item     Var   gamma     se pvalue padj.BH  sig   lower   upper
1 MADRS.S_1 dif.edu -0.0010 0.0246 0.9691  1.0000      -0.0492  0.0473
2 MADRS.S_5 dif.edu -0.0392 0.0228 0.0864  0.5182      -0.0839  0.0056
3 MADRS.S_6 dif.edu -0.1139 0.0238 0.0000  0.0000  *** -0.1606 -0.0672
4 MADRS.S_7 dif.edu  0.0589 0.0256 0.0213  0.1279       0.0088  0.1091
5 MADRS.S_8 dif.edu  0.0770 0.0218 0.0004  0.0025   **  0.0342  0.1198
6 MADRS.S_9 dif.edu  0.0442 0.0256 0.0842  0.5050      -0.0060  0.0944

Code

dif.yr <- dif %>% filter(Time=='WEEK02') %>% select(TreatmentAccessStart) %>% as.vector(.)
dif.yr <- dif.yr[[1]]
RIdifTable(df_r3, dif.yr)

[1] "No statistically significant DIF found."

Code

partgam_DIF(dat.items = as.data.frame(df_r3),
            dat.exo = dif.yr)

       Item    Var   gamma     se pvalue padj.BH  sig   lower   upper
1 MADRS.S_1 dif.yr -0.0037 0.0169 0.8269  1.0000      -0.0367  0.0294
2 MADRS.S_5 dif.yr  0.0177 0.0157 0.2589  1.0000      -0.0130  0.0485
3 MADRS.S_6 dif.yr  0.0030 0.0166 0.8575  1.0000      -0.0296  0.0356
4 MADRS.S_7 dif.yr -0.0105 0.0177 0.5536  1.0000      -0.0452  0.0242
5 MADRS.S_8 dif.yr  0.0297 0.0150 0.0480  0.2880       0.0003  0.0591
6 MADRS.S_9 dif.yr -0.0679 0.0173 0.0001  0.0005  *** -0.1019 -0.0340

11.1 Analysis part 4 decision

Due to the general misfit of infit and outfit (which was the largest) for item 8, aswell as its item-rest-score bootstrap underfitting, removing this.

12 Rasch analysis 5

Removing item 8

Code

df_r3 <- df %>% filter(Time=='WEEK02') %>% select(!Time) %>% select(!c(MADRS.S_4,MADRS.S_3,MADRS.S_2,MADRS.S_8))

Code

RIlistItemsMargin(df_r3, fontsize = 12)

itemnr	item
MADRS.S_1	Sadness
MADRS.S_5	Concentration difficulties
MADRS.S_6	Lassitude
MADRS.S_7	Inability to feel
MADRS.S_9	Suicidal thoughts

Code

RIitemfit(df_r3, cutoff = "Smith98")

Item	InfitMSQ	Location	1 ± 2 / √n
MADRS.S_1	0.889	0.22	[0.973, 1.027]
MADRS.S_5	1.253	-0.35	[0.973, 1.027]
MADRS.S_6	0.943	-0.46	[0.973, 1.027]
MADRS.S_7	0.825	0.08	[0.973, 1.027]
MADRS.S_9	1.113	0.59	[0.973, 1.027]
Note:
MSQ values based on conditional estimation (n = 5465 complete cases).

Code

simfit1 <- RIgetfit(df_r3, iterations = 300, cpu = n_cores)  

RIitemfit(df_r3, simfit1)

Item	InfitMSQ	Infit thresholds	OutfitMSQ	Outfit thresholds	Infit diff	Outfit diff	Location
MADRS.S_1	0.889	[0.953, 1.051]	0.912	[0.952, 1.046]	0.064	0.04	0.22
MADRS.S_5	1.253	[0.96, 1.046]	1.254	[0.962, 1.041]	0.207	0.213	-0.35
MADRS.S_6	0.943	[0.943, 1.044]	0.947	[0.944, 1.043]	no misfit	no misfit	-0.46
MADRS.S_7	0.825	[0.959, 1.048]	0.82	[0.957, 1.046]	0.134	0.137	0.08
MADRS.S_9	1.113	[0.955, 1.039]	1.093	[0.94, 1.057]	0.074	0.036	0.59
Note:
MSQ values based on conditional calculations (n = 5465 complete cases). Simulation based thresholds from 300 simulated datasets.

Code

RIgetfitPlot(simfit1, df_r3)

Code

ICCplot(as.data.frame(df_r3), 
        itemnumber = 1:4, 
        method = "cut", 
        itemdescrip = itemlabels$itemnr[c(1,5,6,7)])

[1] "Please press Zoom on the Plots window to see the plot"

Code

### also suggested:
# library(RASCHplot) # install first with `devtools::install_github("ERRTG/RASCHplot")`
# CICCplot(PCM(df),
#          which.item = c(1:4),
#          lower.groups = c(0,7,14,21,28,35),
#          grid.items = TRUE)

Code

ICCplot(as.data.frame(df_r3), 
        itemnumber = c(5), 
        method = "cut", 
        itemdescrip = itemlabels$itemnr[9])

Code

RIrestscore(df_r3)

Item	Observed value	Model expected value	Absolute difference	Adjusted p-value (BH)	Statistical significance level	Location	Relative location
MADRS.S_1	0.75	0.72	0.03	0.000	***	0.22	2.14
MADRS.S_5	0.66	0.72	0.06	0.000	***	-0.35	1.57
MADRS.S_6	0.75	0.72	0.03	0.000	***	-0.46	1.46
MADRS.S_7	0.78	0.71	0.07	0.000	***	0.08	2.00
MADRS.S_9	0.69	0.71	0.02	0.039	*	0.59	2.51

Code

RIbootRestscore(df_r3,cpu=n_cores,iterations = 500,samplesize=800) # samplesize=600

Item	Item-restscore result	% of iterations	Conditional MSQ infit	Relative average item location
MADRS.S_1	no misfit	41.8	0.89	2.14
MADRS.S_1	overfit	58.2	0.89	2.14
MADRS.S_5	no misfit	12.2	1.25	1.57
MADRS.S_5	underfit	87.8	1.25	1.57
MADRS.S_6	no misfit	73.2	0.94	1.46
MADRS.S_6	overfit	26.8	0.94	1.46
MADRS.S_7	overfit	98.4	0.82	2.00
MADRS.S_9	no misfit	96.4	1.11	2.51
Note:
Results based on 500 bootstrap iterations with a sample size of 800. Conditional mean-square infit based on complete responders only, n = 5465.

Code

RIpcmPCA(df_r3)

PCA of Rasch model residuals
Eigenvalues	Proportion of variance
1.72	36.8%
1.16	23.2%
1.07	21.1%
1.05	18.9%
0.01	0.1%

Code

simcor1 <- RIgetResidCor(df_r3, iterations = 500, cpu = n_cores)
RIresidcorr(df_r3, cutoff = simcor1$p99)

	MADRS.S_1	MADRS.S_5	MADRS.S_6	MADRS.S_7
MADRS.S_1
MADRS.S_5	-0.37
MADRS.S_6	-0.32	-0.06
MADRS.S_7	-0.08	-0.32	-0.2
MADRS.S_9	-0.03	-0.39	-0.34	-0.11
Note:
Relative cut-off value is -0.163, which is 0.059 above the average correlation (-0.222). Correlations above the cut-off are highlighted in red text.

Code

RIpartgamLD(df_r3)

Item 1	Item 2	Partial gamma	SE	Lower CI	Upper CI	Adjusted p-value (BH)
MADRS.S_6	MADRS.S_5	0.324	0.021	0.282	0.366	0.000
MADRS.S_1	MADRS.S_9	0.312	0.023	0.267	0.356	0.000
MADRS.S_7	MADRS.S_9	0.243	0.025	0.194	0.292	0.000
MADRS.S_9	MADRS.S_1	0.242	0.023	0.198	0.287	0.000
MADRS.S_7	MADRS.S_1	0.229	0.024	0.181	0.276	0.000
MADRS.S_5	MADRS.S_6	0.188	0.022	0.145	0.232	0.000
MADRS.S_1	MADRS.S_7	0.176	0.025	0.127	0.225	0.000
MADRS.S_9	MADRS.S_7	0.125	0.025	0.076	0.175	0.000
MADRS.S_7	MADRS.S_6	0.087	0.026	0.037	0.138	0.014

Code

RIloadLoc(df_r3)

Code

mirt(df_r3, model=1, itemtype='Rasch', verbose = FALSE) %>% 
  plot(type="trace", as.table = TRUE, 
       theta_lim = c(-6,6))

Code

# for fewer items or a more magnified figure, use:
#RIitemCats(df)

Code

# increase fig-height above as needed, if you have many items
RItargeting(df_r3)

Code

RIitemHierarchy(df_r3)

Code

iarm::score_groups(as.data.frame(df_r3)) %>% 
  as.data.frame(nm = "score_group") %>% 
  dplyr::count(score_group)

  score_group    n
1           1 2941
2           2 2524

Code

dif_plots <- df_r3 %>% 
  add_column(dif = iarm::score_groups(.)) %>% 
  split(.$dif) %>% # split the data using the DIF variable
  map(~ RItileplot(.x %>% dplyr::select(!dif)) + labs(title = .x$dif))
dif_plots[[1]] + dif_plots[[2]]

Code

clr_tests(df_r3, model = "PCM")


Conditional Likelihood Ratio Tests:

        clr df pvalue sig 
overall 460 28 0       ***

Code

item_obsexp(PCM(df_r3))

Score group 1: 
          mean obs mean exp std.res sig
MADRS.S_1  0.801    0.836   -3.024  -- 
MADRS.S_5  1.293    1.235    4.918  ++ 
MADRS.S_6  1.211    1.216   -0.407     
MADRS.S_7  1.001    1.033   -2.780  -  
MADRS.S_9  0.635    0.622    1.190     

Score group 2: 
          mean obs mean exp std.res sig
MADRS.S_1  2.532    2.496    2.286  +  
MADRS.S_5  2.821    2.882   -3.951  -- 
MADRS.S_6  2.944    2.939    0.316     
MADRS.S_7  2.498    2.465    2.348  +  
MADRS.S_9  1.982    1.996   -1.058

Code

grouping_based_on_score = score_groups(as.data.frame(df_r3))
partgam_DIF(dat.items = as.data.frame(df_r3),
            dat.exo = grouping_based_on_score)

       Item                     Var gamma  se pvalue padj.BH sig lower upper
1 MADRS.S_1 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
2 MADRS.S_5 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
3 MADRS.S_6 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
4 MADRS.S_7 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
5 MADRS.S_9 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA

Code

dif.sex <- dif %>% filter(Time=='WEEK02') %>% select(sex) %>% as.vector(.)
dif.sex <- dif.sex[[1]] #female = 1, male = 2 
RIdifTable(df_r3, dif.sex)

Item	2	3	Mean location	StDev	MaxDiff
MADRS.S_1	0.190	0.232	0.211	0.030	0.042
MADRS.S_5	-0.375	-0.371	-0.373	0.003	0.004
MADRS.S_6	-0.424	-0.579	-0.501	0.110	0.155
MADRS.S_7	0.098	-0.018	0.040	0.082	0.116
MADRS.S_9	0.511	0.736	0.624	0.159	0.224

Code

partgam_DIF(dat.items = as.data.frame(df_r3),
            dat.exo = dif.sex)

       Item     Var   gamma     se pvalue padj.BH sig   lower   upper
1 MADRS.S_1 dif.sex  0.0873 0.0280 0.0019  0.0093  **  0.0323  0.1422
2 MADRS.S_5 dif.sex -0.0571 0.0262 0.0296  0.1479     -0.1085 -0.0057
3 MADRS.S_6 dif.sex -0.0655 0.0281 0.0195  0.0976   . -0.1206 -0.0105
4 MADRS.S_7 dif.sex  0.0987 0.0298 0.0009  0.0046  **  0.0404  0.1571
5 MADRS.S_9 dif.sex -0.0121 0.0287 0.6726  1.0000     -0.0683  0.0441

Code

dif.age <- dif %>% filter(Time=='WEEK02') %>% select(age) %>% as.vector(.)
dif.age <- dif.age[[1]]
RIdifTable(df_r3, dif.age)

Item	3	4	5	Mean location	StDev	MaxDiff
MADRS.S_1	0.208	0.108	0.283	0.200	0.088	0.175
MADRS.S_5	-0.342	-0.296	-0.428	-0.356	0.067	0.132
MADRS.S_6	-0.421	-0.522	-0.470	-0.471	0.050	0.101
MADRS.S_7	0.017	0.145	0.029	0.064	0.070	0.128
MADRS.S_9	0.539	0.566	0.585	0.563	0.023	0.047

Code

partgam_DIF(dat.items = as.data.frame(df_r3),
            dat.exo = dif.age)

       Item     Var   gamma     se pvalue padj.BH  sig   lower   upper
1 MADRS.S_1 dif.age -0.0565 0.0164 0.0006  0.0029   ** -0.0886 -0.0243
2 MADRS.S_5 dif.age -0.0001 0.0153 0.9966  1.0000      -0.0301  0.0299
3 MADRS.S_6 dif.age -0.0556 0.0163 0.0007  0.0034   ** -0.0876 -0.0235
4 MADRS.S_7 dif.age  0.1119 0.0173 0.0000  0.0000  ***  0.0780  0.1458
5 MADRS.S_9 dif.age  0.0491 0.0168 0.0034  0.0172    *  0.0162  0.0819

Code

dif.edu <- dif %>% filter(Time=='WEEK02') %>% select(education) %>% as.vector(.)
dif.edu <- dif.edu[[1]]
RIdifTable(df_r3, dif.edu)

Item	2	3	Mean location	StDev	MaxDiff
MADRS.S_1	0.218	0.176	0.197	0.029	0.041
MADRS.S_5	-0.391	-0.330	-0.360	0.042	0.060
MADRS.S_6	-0.542	-0.365	-0.453	0.125	0.177
MADRS.S_7	0.078	0.042	0.060	0.025	0.036
MADRS.S_9	0.637	0.478	0.557	0.113	0.160

Code

partgam_DIF(dat.items = as.data.frame(df_r3),
            dat.exo = dif.edu)

       Item     Var   gamma     se pvalue padj.BH  sig   lower   upper
1 MADRS.S_1 dif.edu  0.0279 0.0251 0.2671  1.0000      -0.0214  0.0772
2 MADRS.S_5 dif.edu -0.0233 0.0235 0.3210  1.0000      -0.0692  0.0227
3 MADRS.S_6 dif.edu -0.1050 0.0248 0.0000  0.0001  *** -0.1536 -0.0564
4 MADRS.S_7 dif.edu  0.0928 0.0265 0.0005  0.0023   **  0.0408  0.1448
5 MADRS.S_9 dif.edu  0.0589 0.0257 0.0219  0.1096       0.0085  0.1092

Code

dif.yr <- dif %>% filter(Time=='WEEK02') %>% select(TreatmentAccessStart) %>% as.vector(.)
dif.yr <- dif.yr[[1]]
RIdifTable(df_r3, dif.yr)

[1] "No statistically significant DIF found."

Code

partgam_DIF(dat.items = as.data.frame(df_r3),
            dat.exo = dif.yr)

       Item    Var   gamma     se pvalue padj.BH sig   lower   upper
1 MADRS.S_1 dif.yr  0.0144 0.0173 0.4043  1.0000     -0.0195  0.0484
2 MADRS.S_5 dif.yr  0.0318 0.0162 0.0501  0.2504      0.0000  0.0636
3 MADRS.S_6 dif.yr  0.0165 0.0173 0.3402  1.0000     -0.0174  0.0505
4 MADRS.S_7 dif.yr  0.0008 0.0183 0.9653  1.0000     -0.0351  0.0367
5 MADRS.S_9 dif.yr -0.0607 0.0174 0.0005  0.0024  ** -0.0947 -0.0267

12.1 Analysis part 5 decision

5 and 6 seems to load to something else and 5 has the worst fit (infit outfit). Also the partial gamma in relation to 5 and 6. While 7 underfits in the boostrap, the infit/outfit is not as deviant as 5. Removing 5.

13 Rasch analysis 6

Removing item 5

Code

df_r3 <- df %>% filter(Time=='WEEK02') %>% select(!Time) %>% 
  select(!c(MADRS.S_4,MADRS.S_3,MADRS.S_2,MADRS.S_8,MADRS.S_5))

Code

RIlistItemsMargin(df_r3, fontsize = 12)

itemnr	item
MADRS.S_1	Sadness
MADRS.S_6	Lassitude
MADRS.S_7	Inability to feel
MADRS.S_9	Suicidal thoughts

Code

RIitemfit(df_r3, cutoff = "Smith98")

Item	InfitMSQ	Location	1 ± 2 / √n
MADRS.S_1	0.892	0.15	[0.973, 1.027]
MADRS.S_6	1.184	-0.57	[0.973, 1.027]
MADRS.S_7	0.843	-0.01	[0.973, 1.027]
MADRS.S_9	1.11	0.52	[0.973, 1.027]
Note:
MSQ values based on conditional estimation (n = 5465 complete cases).

Code

simfit1 <- RIgetfit(df_r3, iterations = 300, cpu = n_cores)  

RIitemfit(df_r3, simfit1)

Item	InfitMSQ	Infit thresholds	OutfitMSQ	Outfit thresholds	Infit diff	Outfit diff	Location
MADRS.S_1	0.892	[0.953, 1.043]	0.912	[0.95, 1.044]	0.061	0.038	0.15
MADRS.S_6	1.184	[0.959, 1.042]	1.196	[0.962, 1.044]	0.142	0.152	-0.57
MADRS.S_7	0.843	[0.953, 1.057]	0.837	[0.952, 1.054]	0.110	0.115	-0.01
MADRS.S_9	1.11	[0.96, 1.047]	1.096	[0.949, 1.062]	0.063	0.034	0.52
Note:
MSQ values based on conditional calculations (n = 5465 complete cases). Simulation based thresholds from 300 simulated datasets.

Code

RIgetfitPlot(simfit1, df_r3)

Code

ICCplot(as.data.frame(df_r3), 
        itemnumber = 1:4, 
        method = "cut", 
        itemdescrip = itemlabels$itemnr[c(1,6,7,9)])

[1] "Please press Zoom on the Plots window to see the plot"

Code

### also suggested:
# library(RASCHplot) # install first with `devtools::install_github("ERRTG/RASCHplot")`
# CICCplot(PCM(df),
#          which.item = c(1:4),
#          lower.groups = c(0,7,14,21,28,35),
#          grid.items = TRUE)

Code

RIrestscore(df_r3)

Item	Observed value	Model expected value	Absolute difference	Adjusted p-value (BH)	Statistical significance level	Location	Relative location
MADRS.S_1	0.77	0.74	0.03	0.000	***	0.15	2.27
MADRS.S_6	0.71	0.75	0.04	0.000	***	-0.57	1.55
MADRS.S_7	0.80	0.74	0.06	0.000	***	-0.01	2.11
MADRS.S_9	0.72	0.74	0.02	0.042	*	0.52	2.64

Code

RIbootRestscore(df_r3,cpu=n_cores,iterations = 500,samplesize=800) # samplesize=600

Item	Item-restscore result	% of iterations	Conditional MSQ infit	Relative average item location
MADRS.S_1	no misfit	59.8	0.89	2.27
MADRS.S_1	overfit	40.2	0.89	2.27
MADRS.S_6	no misfit	75.0	1.18	1.55
MADRS.S_6	underfit	25.0	1.18	1.55
MADRS.S_7	no misfit	7.8	0.84	2.11
MADRS.S_7	overfit	92.2	0.84	2.11
MADRS.S_9	no misfit	98.0	1.11	2.64
Note:
Results based on 500 bootstrap iterations with a sample size of 800. Conditional mean-square infit based on complete responders only, n = 5465.

Code

RIpcmPCA(df_r3)

PCA of Rasch model residuals
Eigenvalues	Proportion of variance
1.51	41.7%
1.26	30.7%
1.22	27.4%
0.01	0.2%

Code

simcor1 <- RIgetResidCor(df_r3, iterations = 500, cpu = n_cores)
RIresidcorr(df_r3, cutoff = simcor1$p99)

	MADRS.S_1	MADRS.S_6	MADRS.S_7
MADRS.S_1
MADRS.S_6	-0.39
MADRS.S_7	-0.23	-0.25
MADRS.S_9	-0.19	-0.43	-0.26
Note:
Relative cut-off value is -0.225, which is 0.065 above the average correlation (-0.29). Correlations above the cut-off are highlighted in red text.

Code

RIpartgamLD(df_r3)

Item 1	Item 2	Partial gamma	SE	Lower CI	Upper CI	Adjusted p-value (BH)
MADRS.S_7	MADRS.S_6	0.247	0.025	0.198	0.297	0.000
MADRS.S_1	MADRS.S_9	0.216	0.025	0.167	0.265	0.000
MADRS.S_7	MADRS.S_1	0.129	0.027	0.076	0.182	0.000
MADRS.S_7	MADRS.S_9	0.117	0.028	0.063	0.172	0.000
MADRS.S_9	MADRS.S_1	0.092	0.026	0.041	0.142	0.004

Code

RIloadLoc(df_r3)

Code

mirt(df_r3, model=1, itemtype='Rasch', verbose = FALSE) %>% 
  plot(type="trace", as.table = TRUE, 
       theta_lim = c(-6,6))

Code

# for fewer items or a more magnified figure, use:
#RIitemCats(df)

Code

# increase fig-height above as needed, if you have many items
RItargeting(df_r3)

Code

RIitemHierarchy(df_r3)

Code

iarm::score_groups(as.data.frame(df_r3)) %>% 
  as.data.frame(nm = "score_group") %>% 
  dplyr::count(score_group)

  score_group    n
1           1 2965
2           2 2500

Code

dif_plots <- df_r3 %>% 
  add_column(dif = iarm::score_groups(.)) %>% 
  split(.$dif) %>% # split the data using the DIF variable
  map(~ RItileplot(.x %>% dplyr::select(!dif)) + labs(title = .x$dif))
dif_plots[[1]] + dif_plots[[2]]

Code

clr_tests(df_r3, model = "PCM")


Conditional Likelihood Ratio Tests:

        clr df pvalue sig 
overall 316 22 0       ***

Code

item_obsexp(PCM(df_r3))

Score group 1: 
          mean obs mean exp std.res sig
MADRS.S_1  0.812    0.842   -2.778  -  
MADRS.S_6  1.270    1.235    3.255  ++ 
MADRS.S_7  1.026    1.047   -1.968     
MADRS.S_9  0.634    0.618    1.500     

Score group 2: 
          mean obs mean exp std.res sig
MADRS.S_1  2.57     2.54     2.12   +  
MADRS.S_6  2.94     2.98    -2.50   -  
MADRS.S_7  2.52     2.50     1.65      
MADRS.S_9  2.02     2.04    -1.33

Code

grouping_based_on_score = score_groups(as.data.frame(df_r3))
partgam_DIF(dat.items = as.data.frame(df_r3),
            dat.exo = grouping_based_on_score)

       Item                     Var gamma  se pvalue padj.BH sig lower upper
1 MADRS.S_1 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
2 MADRS.S_6 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
3 MADRS.S_7 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
4 MADRS.S_9 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA

Code

dif.sex <- dif %>% filter(Time=='WEEK02') %>% select(sex) %>% as.vector(.)
dif.sex <- dif.sex[[1]] #female = 1, male = 2 
RIdifTable(df_r3, dif.sex)

Item	2	3	Mean location	StDev	MaxDiff
MADRS.S_1	0.111	0.149	0.130	0.027	0.038
MADRS.S_6	-0.550	-0.686	-0.618	0.097	0.136
MADRS.S_7	0.003	-0.112	-0.055	0.082	0.116
MADRS.S_9	0.436	0.650	0.543	0.152	0.214

Code

partgam_DIF(dat.items = as.data.frame(df_r3),
            dat.exo = dif.sex)

       Item     Var   gamma     se pvalue padj.BH sig   lower   upper
1 MADRS.S_1 dif.sex  0.0734 0.0294 0.0125  0.0500   .  0.0158  0.1310
2 MADRS.S_6 dif.sex -0.0837 0.0278 0.0026  0.0105   * -0.1383 -0.0292
3 MADRS.S_7 dif.sex  0.0884 0.0311 0.0044  0.0177   *  0.0275  0.1493
4 MADRS.S_9 dif.sex -0.0391 0.0298 0.1891  0.7564     -0.0974  0.0193

Code

dif.age <- dif %>% filter(Time=='WEEK02') %>% select(age) %>% as.vector(.)
dif.age <- dif.age[[1]]
RIdifTable(df_r3, dif.age)

Item	3	4	5	Mean location	StDev	MaxDiff
MADRS.S_1	0.128	0.049	0.190	0.122	0.071	0.142
MADRS.S_6	-0.527	-0.632	-0.597	-0.585	0.054	0.105
MADRS.S_7	-0.063	0.071	-0.084	-0.025	0.084	0.155
MADRS.S_9	0.462	0.513	0.490	0.489	0.025	0.051

Code

partgam_DIF(dat.items = as.data.frame(df_r3),
            dat.exo = dif.age)

       Item     Var   gamma     se pvalue padj.BH  sig   lower   upper
1 MADRS.S_1 dif.age -0.0688 0.0172 0.0001  0.0003  *** -0.1025 -0.0351
2 MADRS.S_6 dif.age -0.0539 0.0162 0.0009  0.0035   ** -0.0856 -0.0222
3 MADRS.S_7 dif.age  0.1161 0.0179 0.0000  0.0000  ***  0.0810  0.1511
4 MADRS.S_9 dif.age  0.0522 0.0175 0.0029  0.0117    *  0.0178  0.0866

Code

dif.edu <- dif %>% filter(Time=='WEEK02') %>% select(education) %>% as.vector(.)
dif.edu <- dif.edu[[1]]
RIdifTable(df_r3, dif.edu)

Item	2	3	Mean location	StDev	MaxDiff
MADRS.S_1	0.133	0.109	0.121	0.017	0.023
MADRS.S_6	-0.670	-0.469	-0.570	0.142	0.201
MADRS.S_7	-0.019	-0.050	-0.034	0.022	0.032
MADRS.S_9	0.556	0.410	0.483	0.103	0.146

Code

partgam_DIF(dat.items = as.data.frame(df_r3),
            dat.exo = dif.edu)

       Item     Var   gamma     se pvalue padj.BH  sig   lower   upper
1 MADRS.S_1 dif.edu  0.0217 0.0263 0.4104  1.0000      -0.0299  0.0732
2 MADRS.S_6 dif.edu -0.1160 0.0245 0.0000  0.0000  *** -0.1640 -0.0679
3 MADRS.S_7 dif.edu  0.0811 0.0276 0.0033  0.0132    *  0.0270  0.1351
4 MADRS.S_9 dif.edu  0.0713 0.0267 0.0076  0.0304    *  0.0190  0.1237

Code

dif.yr <- dif %>% filter(Time=='WEEK02') %>% select(TreatmentAccessStart) %>% as.vector(.)
dif.yr <- dif.yr[[1]]
RIdifTable(df_r3, dif.yr)

[1] "No statistically significant DIF found."

Code

partgam_DIF(dat.items = as.data.frame(df_r3),
            dat.exo = dif.yr)

       Item    Var   gamma     se pvalue padj.BH sig   lower   upper
1 MADRS.S_1 dif.yr  0.0235 0.0179 0.1899  0.7595     -0.0116  0.0587
2 MADRS.S_6 dif.yr  0.0293 0.0172 0.0887  0.3548     -0.0044  0.0631
3 MADRS.S_7 dif.yr  0.0063 0.0189 0.7403  1.0000     -0.0308  0.0434
4 MADRS.S_9 dif.yr -0.0611 0.0181 0.0007  0.0029  ** -0.0966 -0.0257

13.1 Analysis part 6 decision

7 or 6, trying to remove both in separate analyses.

14 Rasch analysis 7 part 1 of 2

Removing item 6

Code

df_r3 <- df %>% filter(Time=='WEEK02') %>% select(!Time) %>% 
  select(!c(MADRS.S_4,MADRS.S_3,MADRS.S_2,MADRS.S_8,MADRS.S_5,MADRS.S_6))

Code

RIlistItemsMargin(df_r3, fontsize = 12)

itemnr	item
MADRS.S_1	Sadness
MADRS.S_7	Inability to feel
MADRS.S_9	Suicidal thoughts

Code

RIitemfit(df_r3, cutoff = "Smith98")

Item	InfitMSQ	Location	1 ± 2 / √n
MADRS.S_1	0.943	-0.05	[0.973, 1.027]
MADRS.S_7	0.971	-0.22	[0.973, 1.027]
MADRS.S_9	1.096	0.33	[0.973, 1.027]
Note:
MSQ values based on conditional estimation (n = 5465 complete cases).

Code

simfit1 <- RIgetfit(df_r3, iterations = 300, cpu = n_cores)  

RIitemfit(df_r3, simfit1)

Item	InfitMSQ	Infit thresholds	OutfitMSQ	Outfit thresholds	Infit diff	Outfit diff	Location
MADRS.S_1	0.943	[0.96, 1.039]	0.956	[0.96, 1.04]	0.017	0.004	-0.05
MADRS.S_7	0.971	[0.967, 1.033]	0.967	[0.966, 1.032]	no misfit	no misfit	-0.22
MADRS.S_9	1.096	[0.956, 1.033]	1.086	[0.954, 1.038]	0.063	0.048	0.33
Note:
MSQ values based on conditional calculations (n = 5465 complete cases). Simulation based thresholds from 300 simulated datasets.

Code

RIgetfitPlot(simfit1, df_r3)

Code

ICCplot(as.data.frame(df_r3), 
        itemnumber = 1:3, 
        method = "cut", 
        itemdescrip = itemlabels$itemnr[c(1,7,9)])

[1] "Please press Zoom on the Plots window to see the plot"

Code

### also suggested:
# library(RASCHplot) # install first with `devtools::install_github("ERRTG/RASCHplot")`
# CICCplot(PCM(df),
#          which.item = c(1:4),
#          lower.groups = c(0,7,14,21,28,35),
#          grid.items = TRUE)

Code

RIrestscore(df_r3)

Item	Observed value	Model expected value	Absolute difference	Adjusted p-value (BH)	Statistical significance level	Location	Relative location
MADRS.S_1	0.78	0.76	0.02	0.023	*	-0.05	2.39
MADRS.S_7	0.78	0.76	0.02	0.044	*	-0.22	2.22
MADRS.S_9	0.75	0.76	0.01	0.119		0.33	2.77

Code

RIbootRestscore(df_r3,cpu=n_cores,iterations = 500,samplesize=800) # samplesize=600

Item	Item-restscore result	% of iterations	Conditional MSQ infit	Relative average item location
MADRS.S_1	no misfit	96.2	0.94	2.39
MADRS.S_7	no misfit	98.4	0.97	2.22
MADRS.S_9	no misfit	100.0	1.10	2.77
Note:
Results based on 500 bootstrap iterations with a sample size of 800. Conditional mean-square infit based on complete responders only, n = 5465.

Code

RIpcmPCA(df_r3)

PCA of Rasch model residuals
Eigenvalues	Proportion of variance
1.51	53.1%
1.49	46.6%
0.01	0.3%

Code

simcor1 <- RIgetResidCor(df_r3, iterations = 500, cpu = n_cores)
RIresidcorr(df_r3, cutoff = simcor1$p99)

	MADRS.S_1	MADRS.S_7
MADRS.S_1
MADRS.S_7	-0.37
MADRS.S_9	-0.43	-0.43
Note:
Relative cut-off value is -0.354, which is 0.06 above the average correlation (-0.413). Correlations above the cut-off are highlighted in red text.

Code

RIpartgamLD(df_r3)

Item 1	Item 2	Partial gamma	SE	Lower CI	Upper CI	Adjusted p-value (BH)
MADRS.S_1	MADRS.S_7	0.11	0.03	0.051	0.170	0.002
MADRS.S_7	MADRS.S_1	0.103	0.03	0.044	0.162	0.004

Code

RIloadLoc(df_r3)

Code

mirt(df_r3, model=1, itemtype='Rasch', verbose = FALSE) %>% 
  plot(type="trace", as.table = TRUE, 
       theta_lim = c(-6,6))

Code

# for fewer items or a more magnified figure, use:
#RIitemCats(df)

Code

# increase fig-height above as needed, if you have many items
RItargeting(df_r3)

Code

RIitemHierarchy(df_r3)

Code

iarm::score_groups(as.data.frame(df_r3)) %>% 
  as.data.frame(nm = "score_group") %>% 
  dplyr::count(score_group)

  score_group    n
1           1 2961
2           2 2504

Code

dif_plots <- df_r3 %>% 
  add_column(dif = iarm::score_groups(.)) %>% 
  split(.$dif) %>% # split the data using the DIF variable
  map(~ RItileplot(.x %>% dplyr::select(!dif)) + labs(title = .x$dif))
dif_plots[[1]] + dif_plots[[2]]

Code

clr_tests(df_r3, model = "PCM")


Conditional Likelihood Ratio Tests:

        clr df pvalue sig 
overall 134 16 0       ***

Code

item_obsexp(PCM(df_r3))

Score group 1: 
          mean obs mean exp std.res sig
MADRS.S_1  0.868    0.894   -2.438  -  
MADRS.S_7  1.129    1.124    0.507     
MADRS.S_9  0.661    0.641    1.924     

Score group 2: 
          mean obs mean exp std.res sig
MADRS.S_1  2.592    2.568    1.927     
MADRS.S_7  2.526    2.531   -0.423     
MADRS.S_9  2.054    2.074   -1.701

Code

grouping_based_on_score = score_groups(as.data.frame(df_r3))
partgam_DIF(dat.items = as.data.frame(df_r3),
            dat.exo = grouping_based_on_score)

       Item                     Var gamma  se pvalue padj.BH sig lower upper
1 MADRS.S_1 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
2 MADRS.S_7 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
3 MADRS.S_9 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA

Code

dif.sex <- dif %>% filter(Time=='WEEK02') %>% select(sex) %>% as.vector(.)
dif.sex <- dif.sex[[1]] #female = 1, male = 2 
RIdifTable(df_r3, dif.sex)

[1] "No statistically significant DIF found."

Code

partgam_DIF(dat.items = as.data.frame(df_r3),
            dat.exo = dif.sex)

       Item     Var   gamma     se pvalue padj.BH sig   lower   upper
1 MADRS.S_1 dif.sex  0.0377 0.0315 0.2313  0.6939     -0.0240  0.0995
2 MADRS.S_7 dif.sex  0.0597 0.0322 0.0637  0.1911     -0.0034  0.1229
3 MADRS.S_9 dif.sex -0.0840 0.0319 0.0084  0.0252   * -0.1465 -0.0215

Code

dif.age <- dif %>% filter(Time=='WEEK02') %>% select(age) %>% as.vector(.)
dif.age <- dif.age[[1]]
RIdifTable(df_r3, dif.age)

Item	3	4	5	Mean location	StDev	MaxDiff
MADRS.S_1	-0.058	-0.169	0.007	-0.073	0.089	0.176
MADRS.S_7	-0.243	-0.148	-0.304	-0.232	0.079	0.156
MADRS.S_9	0.301	0.317	0.297	0.305	0.011	0.020

Code

partgam_DIF(dat.items = as.data.frame(df_r3),
            dat.exo = dif.age)

       Item     Var   gamma     se pvalue padj.BH  sig   lower   upper
1 MADRS.S_1 dif.age -0.1192 0.0181 0.0000  0.0000  *** -0.1546 -0.0838
2 MADRS.S_7 dif.age  0.0986 0.0184 0.0000  0.0000  ***  0.0625  0.1348
3 MADRS.S_9 dif.age  0.0290 0.0189 0.1259  0.3778      -0.0081  0.0661

Code

dif.edu <- dif %>% filter(Time=='WEEK02') %>% select(education) %>% as.vector(.)
dif.edu <- dif.edu[[1]]
RIdifTable(df_r3, dif.edu)

Item	2	3	Mean location	StDev	MaxDiff
MADRS.S_1	-0.086	-0.042	-0.064	0.031	0.044
MADRS.S_7	-0.254	-0.223	-0.239	0.022	0.031
MADRS.S_9	0.340	0.265	0.303	0.054	0.076

Code

partgam_DIF(dat.items = as.data.frame(df_r3),
            dat.exo = dif.edu)

       Item     Var   gamma     se pvalue padj.BH sig   lower  upper
1 MADRS.S_1 dif.edu -0.0443 0.0280 0.1137  0.3412     -0.0993 0.0106
2 MADRS.S_7 dif.edu  0.0240 0.0287 0.4028  1.0000     -0.0322 0.0803
3 MADRS.S_9 dif.edu  0.0120 0.0288 0.6775  1.0000     -0.0444 0.0683

Code

dif.yr <- dif %>% filter(Time=='WEEK02') %>% select(TreatmentAccessStart) %>% as.vector(.)
dif.yr <- dif.yr[[1]]
RIdifTable(df_r3, dif.yr)

[1] "No statistically significant DIF found."

Code

partgam_DIF(dat.items = as.data.frame(df_r3),
            dat.exo = dif.yr)

       Item    Var   gamma     se pvalue padj.BH sig   lower   upper
1 MADRS.S_1 dif.yr  0.0438 0.0192 0.0222  0.0667   .  0.0063  0.0814
2 MADRS.S_7 dif.yr  0.0233 0.0195 0.2328  0.6985     -0.0150  0.0615
3 MADRS.S_9 dif.yr -0.0674 0.0194 0.0005  0.0015  ** -0.1054 -0.0294

14.1 Analysis part 7 1 of 2

This was removing item 6

15 Rasch analysis 7 part 2 of 2

Removing item 7

Code

df_r3 <- df %>% filter(Time=='WEEK02') %>% select(!Time) %>% 
  select(!c(MADRS.S_4,MADRS.S_3,MADRS.S_2,MADRS.S_8,MADRS.S_5,MADRS.S_7))

Code

RIlistItemsMargin(df_r3, fontsize = 12)

itemnr	item
MADRS.S_1	Sadness
MADRS.S_6	Lassitude
MADRS.S_9	Suicidal thoughts

Code

RIitemfit(df_r3, cutoff = "Smith98")

Item	InfitMSQ	Location	1 ± 2 / √n
MADRS.S_1	0.857	0.13	[0.973, 1.027]
MADRS.S_6	1.137	-0.54	[0.973, 1.027]
MADRS.S_9	1.031	0.50	[0.973, 1.027]
Note:
MSQ values based on conditional estimation (n = 5465 complete cases).

Code

simfit1 <- RIgetfit(df_r3, iterations = 300, cpu = n_cores)  

RIitemfit(df_r3, simfit1)

Item	InfitMSQ	Infit thresholds	OutfitMSQ	Outfit thresholds	Infit diff	Outfit diff	Location
MADRS.S_1	0.857	[0.966, 1.032]	0.869	[0.962, 1.033]	0.109	0.093	0.13
MADRS.S_6	1.137	[0.958, 1.034]	1.131	[0.961, 1.036]	0.103	0.095	-0.54
MADRS.S_9	1.031	[0.96, 1.038]	1.017	[0.945, 1.038]	no misfit	no misfit	0.50
Note:
MSQ values based on conditional calculations (n = 5465 complete cases). Simulation based thresholds from 300 simulated datasets.

Code

RIgetfitPlot(simfit1, df_r3)

Code

ICCplot(as.data.frame(df_r3), 
        itemnumber = 1:3, 
        method = "cut", 
        itemdescrip = itemlabels$itemnr[c(1,6,9)])

[1] "Please press Zoom on the Plots window to see the plot"

Code

### also suggested:
# library(RASCHplot) # install first with `devtools::install_github("ERRTG/RASCHplot")`
# CICCplot(PCM(df),
#          which.item = c(1:4),
#          lower.groups = c(0,7,14,21,28,35),
#          grid.items = TRUE)

Code

RIrestscore(df_r3)

Item	Observed value	Model expected value	Absolute difference	Adjusted p-value (BH)	Statistical significance level	Location	Relative location
MADRS.S_1	0.76	0.71	0.05	0.000	***	0.13	2.10
MADRS.S_6	0.68	0.71	0.03	0.001	**	-0.54	1.42
MADRS.S_9	0.71	0.70	0.01	0.466		0.50	2.46

Code

RIbootRestscore(df_r3,cpu=n_cores,iterations = 500,samplesize=800) # samplesize=600

Item	Item-restscore result	% of iterations	Conditional MSQ infit	Relative average item location
MADRS.S_1	no misfit	32.0	0.86	2.10
MADRS.S_1	overfit	68.0	0.86	2.10
MADRS.S_6	no misfit	97.6	1.14	1.42
MADRS.S_9	no misfit	100.0	1.03	2.46
Note:
Results based on 500 bootstrap iterations with a sample size of 800. Conditional mean-square infit based on complete responders only, n = 5465.

Code

RIpcmPCA(df_r3)

PCA of Rasch model residuals
Eigenvalues	Proportion of variance
1.66	57.6%
1.34	42.1%
0.01	0.2%

Code

simcor1 <- RIgetResidCor(df_r3, iterations = 500, cpu = n_cores)
RIresidcorr(df_r3, cutoff = simcor1$p99)

	MADRS.S_1	MADRS.S_6
MADRS.S_1
MADRS.S_6	-0.45
MADRS.S_9	-0.26	-0.51
Note:
Relative cut-off value is -0.348, which is 0.06 above the average correlation (-0.408). Correlations above the cut-off are highlighted in red text.

Code

RIpartgamLD(df_r3)

Item 1	Item 2	Partial gamma	SE	Lower CI	Upper CI
MADRS.S_1	MADRS.S_9	0.275	0.026	0.224	0.326
MADRS.S_1	MADRS.S_6	0.16	0.028	0.106	0.214
MADRS.S_9	MADRS.S_1	0.116	0.027	0.064	0.169

Code

RIloadLoc(df_r3)

Code

mirt(df_r3, model=1, itemtype='Rasch', verbose = FALSE) %>% 
  plot(type="trace", as.table = TRUE, 
       theta_lim = c(-6,6))

Code

# for fewer items or a more magnified figure, use:
#RIitemCats(df)

Code

# increase fig-height above as needed, if you have many items
RItargeting(df_r3)

Code

RIitemHierarchy(df_r3)

Code

iarm::score_groups(as.data.frame(df_r3)) %>% 
  as.data.frame(nm = "score_group") %>% 
  dplyr::count(score_group)

  score_group    n
1           1 2816
2           2 2649

Code

dif_plots <- df_r3 %>% 
  add_column(dif = iarm::score_groups(.)) %>% 
  split(.$dif) %>% # split the data using the DIF variable
  map(~ RItileplot(.x %>% dplyr::select(!dif)) + labs(title = .x$dif))
dif_plots[[1]] + dif_plots[[2]]

Code

clr_tests(df_r3, model = "PCM")


Conditional Likelihood Ratio Tests:

        clr df pvalue sig 
overall 243 16 0       ***

Code

item_obsexp(PCM(df_r3))

Score group 1: 
          mean obs mean exp std.res sig
MADRS.S_1  0.780    0.819   -3.489  -- 
MADRS.S_6  1.252    1.222    2.763  +  
MADRS.S_9  0.593    0.585    0.735     

Score group 2: 
          mean obs mean exp std.res sig
MADRS.S_1  2.515    2.480    2.509  +  
MADRS.S_6  2.894    2.921   -1.975     
MADRS.S_9  1.992    2.000   -0.595

Code

grouping_based_on_score = score_groups(as.data.frame(df_r3))
partgam_DIF(dat.items = as.data.frame(df_r3),
            dat.exo = grouping_based_on_score)

       Item                     Var gamma  se pvalue padj.BH sig lower upper
1 MADRS.S_1 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
2 MADRS.S_6 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA
3 MADRS.S_9 grouping_based_on_score   NaN NaN     NA      NA   ?    NA    NA

Code

dif.sex <- dif %>% filter(Time=='WEEK02') %>% select(sex) %>% as.vector(.)
dif.sex <- dif.sex[[1]] #female = 1, male = 2 
RIdifTable(df_r3, dif.sex)

Item	2	3	Mean location	StDev	MaxDiff
MADRS.S_1	0.101	0.102	0.101	0.001	0.001
MADRS.S_6	-0.513	-0.701	-0.607	0.133	0.188
MADRS.S_9	0.412	0.599	0.506	0.133	0.188

Code

partgam_DIF(dat.items = as.data.frame(df_r3),
            dat.exo = dif.sex)

       Item     Var   gamma     se pvalue padj.BH  sig   lower   upper
1 MADRS.S_1 dif.sex  0.1089 0.0301 0.0003  0.0009  ***  0.0499  0.1680
2 MADRS.S_6 dif.sex -0.0674 0.0286 0.0186  0.0557    . -0.1234 -0.0113
3 MADRS.S_9 dif.sex -0.0253 0.0307 0.4085  1.0000      -0.0854  0.0347

Code

dif.age <- dif %>% filter(Time=='WEEK02') %>% select(age) %>% as.vector(.)
dif.age <- dif.age[[1]]
RIdifTable(df_r3, dif.age)

Item	2	3	Mean location	StDev	MaxDiff
MADRS.S_1	0.067	0.121	0.094	0.039	0.055
MADRS.S_6	-0.626	-0.554	-0.590	0.051	0.072
MADRS.S_9	0.559	0.433	0.496	0.090	0.127

Code

partgam_DIF(dat.items = as.data.frame(df_r3),
            dat.exo = dif.age)

       Item     Var   gamma     se pvalue padj.BH  sig   lower   upper
1 MADRS.S_1 dif.age -0.0434 0.0177 0.0143  0.0430    * -0.0782 -0.0087
2 MADRS.S_6 dif.age -0.0271 0.0166 0.1021  0.3062      -0.0596  0.0054
3 MADRS.S_9 dif.age  0.0885 0.0179 0.0000  0.0000  ***  0.0534  0.1236

Code

dif.edu <- dif %>% filter(Time=='WEEK02') %>% select(education) %>% as.vector(.)
dif.edu <- dif.edu[[1]]
RIdifTable(df_r3, dif.edu)

Item	2	3	Mean location	StDev	MaxDiff
MADRS.S_1	0.112	0.083	0.098	0.021	0.030
MADRS.S_6	-0.642	-0.454	-0.548	0.133	0.189
MADRS.S_9	0.530	0.371	0.450	0.112	0.159

Code

partgam_DIF(dat.items = as.data.frame(df_r3),
            dat.exo = dif.edu)

       Item     Var   gamma     se pvalue padj.BH  sig   lower   upper
1 MADRS.S_1 dif.edu  0.0468 0.0270 0.0828  0.2484      -0.0061  0.0998
2 MADRS.S_6 dif.edu -0.1020 0.0253 0.0001  0.0002  *** -0.1515 -0.0524
3 MADRS.S_9 dif.edu  0.0981 0.0272 0.0003  0.0009  ***  0.0448  0.1514

Code

dif.yr <- dif %>% filter(Time=='WEEK02') %>% select(TreatmentAccessStart) %>% as.vector(.)
dif.yr <- dif.yr[[1]]
RIdifTable(df_r3, dif.yr)

Item	3	4	5	Mean location	StDev	MaxDiff
MADRS.S_1	0.207	-0.387	0.043	-0.046	0.307	0.594
MADRS.S_6	-0.643	-0.068	-0.568	-0.426	0.312	0.574
MADRS.S_9	0.436	0.455	0.525	0.472	0.047	0.089

Code

partgam_DIF(dat.items = as.data.frame(df_r3),
            dat.exo = dif.yr)

       Item    Var   gamma     se pvalue padj.BH  sig   lower   upper
1 MADRS.S_1 dif.yr  0.0309 0.0183 0.0918  0.2755      -0.0050  0.0669
2 MADRS.S_6 dif.yr  0.0341 0.0175 0.0517  0.1551      -0.0002  0.0685
3 MADRS.S_9 dif.yr -0.0688 0.0186 0.0002  0.0006  *** -0.1052 -0.0324

15.1 Analysis part 7 2 of 2

This was removing item 7

15.2 Analysis part 7 conclusion

Removing item 6 resulted in a better fitting scale.

16 Analysis conclusion

Resulting MADRS-S reworked scale thus has item 1,7 and 9.

16.1 Resulting items

Code

df_r3 <- df %>% filter(Time=='WEEK02') %>% select(!Time) %>% 
  select(!c(MADRS.S_4,MADRS.S_3,MADRS.S_2,MADRS.S_8,MADRS.S_5,MADRS.S_6))

what_scale <- gsub('\\.','',items_to_use)

Code

RIlistItemsMargin(df_r3, fontsize = 12)

itemnr	item
MADRS.S_1	Sadness
MADRS.S_7	Inability to feel
MADRS.S_9	Suicidal thoughts

16.2 Reliability

Code

RItif(df_r3)

Code

RItif(df_r3,samplePSI=TRUE)

16.3 Person parameter

Code

RIpfit(df_r3)

16.4 Misfit check

Code

pfit_u3poly <- PerFit::U3poly(matrix = df_r3, 
                      Ncat = 7, # make sure to input number of response categories, not thresholds
                      IRT.PModel = "PCM")
misfits <- PerFit::flagged.resp(pfit_u3poly) %>% 
  pluck("Scores") %>% 
  as.data.frame() %>% 
  pull(FlaggedID)

misfits2 <- RIpfit(df_r3, output = "rowid")

Code

RIitemfit(df_r3, simcut = simfit1)

Item	InfitMSQ	Infit thresholds	OutfitMSQ	Outfit thresholds	Infit diff	Outfit diff	Location
MADRS.S_1	0.943	[0.966, 1.032]	0.956	[0.962, 1.033]	0.023	0.006	-0.05
MADRS.S_7	0.971	[0.958, 1.034]	0.967	[0.961, 1.036]	no misfit	no misfit	-0.22
MADRS.S_9	1.096	[0.96, 1.038]	1.086	[0.945, 1.038]	0.058	0.048	0.33
Note:
MSQ values based on conditional calculations (n = 5465 complete cases). Simulation based thresholds from 300 simulated datasets.

Code

RIitemfit(df_r3[-misfits,], simcut = simfit1)

Item	InfitMSQ	Infit thresholds	OutfitMSQ	Outfit thresholds	Infit diff	Outfit diff	Location
MADRS.S_1	0.953	[0.966, 1.032]	0.904	[0.962, 1.033]	0.013	0.058	0.12
MADRS.S_7	0.984	[0.958, 1.034]	0.862	[0.961, 1.036]	no misfit	0.099	-0.42
MADRS.S_9	1.085	[0.96, 1.038]	1.028	[0.945, 1.038]	0.047	no misfit	0.38
Note:
MSQ values based on conditional calculations (n = 5006 complete cases). Simulation based thresholds from 300 simulated datasets.

Code

RItif(df_r3[-misfits,])

Code

RIitemfit(df_r3[-misfits2,], simcut = simfit1)

Item	InfitMSQ	Infit thresholds	OutfitMSQ	Outfit thresholds	Infit diff	Outfit diff	Location
MADRS.S_1	0.958	[0.966, 1.032]	0.96	[0.962, 1.033]	0.008	0.002	-0.09
MADRS.S_7	0.982	[0.958, 1.034]	0.973	[0.961, 1.036]	no misfit	no misfit	-0.20
MADRS.S_9	1.064	[0.96, 1.038]	1.047	[0.945, 1.038]	0.026	0.009	0.35
Note:
MSQ values based on conditional calculations (n = 5229 complete cases). Simulation based thresholds from 300 simulated datasets.

Code

RItif(df_r3[-misfits2,])

16.5 Item parameters

Code

RIitemparams(df_r3)

	Threshold 1	Threshold 2	Threshold 3	Threshold 4	Threshold 5	Threshold 6	Item location
MADRS.S_1	-4.46	-2.08	-0.01	0.38	2.35	3.52	-0.05
MADRS.S_7	-5.08	-2.60	-0.14	1.28	2.38	2.82	-0.22
MADRS.S_9	-3.54	-1.93	0.96	2.11	4.03	NA	0.33
Note:
Item location is the average of the thresholds for each item.

Code

item_param_matrix <- RIitemparams(df_r3,output='dataframe') %>% as.matrix(.)
write.csv(item_param_matrix,file=paste0("./results/item_params_",what_scale,'.csv'),row.names=TRUE)

16.6 Ordinal sum to interval score

Code

RIscoreSE(df_r3)

Ordinal sum score	Logit score	Logit std.error
0	-6.470	0.789
1	-4.981	1.002
2	-4.021	0.951
3	-3.219	0.893
4	-2.494	0.867
5	-1.767	0.846
6	-1.012	0.811
7	-0.358	0.771
8	0.174	0.735
9	0.659	0.703
10	1.131	0.680
11	1.586	0.669
12	2.016	0.672
13	2.437	0.691
14	2.880	0.735
15	3.401	0.805
16	4.112	0.869
17	5.444	0.719

Code

RIscoreSE(df_r3,output='figure')

Code

sum_to_latent <- RIscoreSE(df_r3,output = 'dataframe')
write.csv(sum_to_latent,file=paste0('./results/ordinal_sum_to_latent_',what_scale,'.csv'))

16.7 Thetas

The ordinal sum to interval score contains the information to transform into the below thetas, but we plot them here for convinience (based on 3 items)

Code

est_thetas2 <- RIestThetas(df_r3, method = "WLE")
hist(est_thetas2$WLE, 
     col = "#009ca6", 
     main = "Histogram of person locations (thetas)", 
     breaks = 17)

17 Software used

Code

pkgs <- grateful::cite_packages(cite.tidyverse = TRUE, 
                      output = "table",
                      bib.file = "grateful-refs.bib",
                      include.RStudio = FALSE,
                      out.dir = getwd())
formattable(pkgs, table.attr = 'class=\"table table-striped\" style="font-size: 15px; font-family: Lato; width: 80%"')

Package	Version	Citation
base	4.4.1	@base
car	3.1.3	@car
doParallel	1.0.17	@doParallel
easyRasch	0.3.3	@easyRasch
eRm	1.0.6	@eRm2007a; @eRm2007b; @eRm2009c; @eRm2013d; @eRm2015e; @eRm2019f
formattable	0.2.1	@formattable
furrr	0.3.1	@furrr
ggrepel	0.9.6	@ggrepel
glue	1.8.0	@glue
gridExtra	2.3	@gridExtra
iarm	0.4.3	@iarm
janitor	2.2.0	@janitor
kableExtra	1.4.0	@kableExtra
knitr	1.49	@knitr2014; @knitr2015; @knitr2024
matrixStats	1.4.1	@matrixStats
mirt	1.43	@mirt
mokken	3.1.2	@mokken2007; @mokken2012
patchwork	1.3.0	@patchwork
PerFit	1.4.6	@PerFit
psych	2.4.6.26	@psych
psychotools	0.7.4	@psychotools2021; @psychotools2022; @psychotools2024
psychotree	0.16.1	@psychotree2010e; @psychotree2011a; @psychotree2015b; @psychotree2018c; @psychotree2018d
reshape	0.8.9	@reshape
rmarkdown	2.29	@rmarkdown2018; @rmarkdown2020; @rmarkdown2024
tidyverse	2.0.0	@tidyverse

18 Session info

Code

sessionInfo()

R version 4.4.1 (2024-06-14)
Platform: x86_64-pc-linux-gnu
Running under: Ubuntu 22.04.5 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.10.0 
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

time zone: Europe/Stockholm
tzcode source: system (glibc)

attached base packages:
 [1] parallel  grid      stats4    stats     graphics  grDevices utils    
 [8] datasets  methods   base     

other attached packages:
 [1] gridExtra_2.3     knitr_1.49        car_3.1-3         carData_3.0-5    
 [5] grateful_0.2.10   easyRasch_0.3.3   doParallel_1.0.17 iterators_1.0.14 
 [9] furrr_0.3.1       future_1.34.0     foreach_1.5.2     ggdist_3.3.2     
[13] janitor_2.2.0     iarm_0.4.3        hexbin_1.28.5     catR_3.17        
[17] glue_1.8.0        ggrepel_0.9.6     patchwork_1.3.0   reshape_0.8.9    
[21] matrixStats_1.4.1 psychotree_0.16-1 psychotools_0.7-4 partykit_1.2-23  
[25] mvtnorm_1.3-2     libcoin_1.0-10    psych_2.4.6.26    mirt_1.43        
[29] lattice_0.22-5    eRm_1.0-6         lubridate_1.9.4   forcats_1.0.0    
[33] stringr_1.5.1     dplyr_1.1.4       purrr_1.0.2       readr_2.1.5      
[37] tidyr_1.3.1       tibble_3.2.1      ggplot2_3.5.1     tidyverse_2.0.0  
[41] kableExtra_1.4.0  formattable_0.2.1

loaded via a namespace (and not attached):
  [1] splines_4.4.1        bitops_1.0-6         R.oo_1.27.0         
  [4] cellranger_1.1.0     deSolve_1.40         rpart_4.1.23        
  [7] lifecycle_1.0.4      rprojroot_2.0.4      globals_0.16.3      
 [10] MASS_7.3-61          backports_1.5.0      magrittr_2.0.3      
 [13] vcd_1.4-13           Hmisc_5.2-1          rmarkdown_2.29      
 [16] yaml_2.3.10          sm_2.2-6.0           sessioninfo_1.2.2   
 [19] pbapply_1.7-2        RColorBrewer_1.1-3   abind_1.4-8         
 [22] audio_0.1-11         expm_1.0-0           quadprog_1.5-8      
 [25] R.utils_2.12.3       RCurl_1.98-1.16      pracma_2.4.4        
 [28] nnet_7.3-19          listenv_0.9.1        testthat_3.2.2      
 [31] RPushbullet_0.3.4    vegan_2.6-8          parallelly_1.40.1   
 [34] svglite_2.1.3        permute_0.9-7        codetools_0.2-19    
 [37] xml2_1.3.6           tidyselect_1.2.1     farver_2.1.2        
 [40] base64enc_0.1-3      jsonlite_1.8.9       polycor_0.8-1       
 [43] ks_1.14.3            progressr_0.15.1     Formula_1.2-5       
 [46] survival_3.7-0       systemfonts_1.1.0    tools_4.4.1         
 [49] gnm_1.1-5            ragg_1.3.3           snow_0.4-4          
 [52] Rcpp_1.0.13-1        mnormt_2.1.1         admisc_0.36         
 [55] xfun_0.49            here_1.0.1           mgcv_1.9-1          
 [58] msm_1.8.2            distributional_0.5.0 ca_0.71.1           
 [61] withr_3.0.2          beepr_2.0            fastmap_1.2.0       
 [64] fansi_1.0.6          digest_0.6.37        timechange_0.3.0    
 [67] R6_2.5.1             textshaping_0.4.1    colorspace_2.1-1    
 [70] R.methodsS3_1.8.2    inum_1.0-5           utf8_1.2.4          
 [73] generics_0.1.3       renv_1.0.10          data.table_1.16.4   
 [76] SimDesign_2.17.1     htmlwidgets_1.6.4    pkgconfig_2.0.3     
 [79] gtable_0.3.6         lmtest_0.9-40        pcaPP_2.0-5         
 [82] brio_1.1.5           htmltools_0.5.8.1    scales_1.3.0        
 [85] snakecase_0.11.1     irtoys_0.2.2         rstudioapi_0.17.1   
 [88] tzdb_0.4.0           checkmate_2.3.2      nlme_3.1-165        
 [91] curl_6.0.1           zoo_1.8-12           KernSmooth_2.23-24  
 [94] relimp_1.0-5         vcdExtra_0.8-5       fda_6.2.0           
 [97] foreign_0.8-86       pillar_1.9.0         vctrs_0.6.5         
[100] Deriv_4.1.6          cluster_2.1.6        dcurver_0.9.2       
[103] GPArotation_2024.3-1 htmlTable_2.4.3      evaluate_1.0.1      
[106] cli_3.6.3            compiler_4.4.1       rlang_1.1.4         
[109] future.apply_1.11.3  labeling_0.4.3       mclust_6.1.1        
[112] fds_1.8              plyr_1.8.9           rainbow_3.8         
[115] stringi_1.8.4        hdrcde_3.4           viridisLite_0.4.2   
[118] munsell_0.5.1        Matrix_1.6-5         qvcalc_1.0.3        
[121] hms_1.1.3            ltm_1.2-0            PerFit_1.4.6        
[124] readxl_1.4.3

19 References

Reuse

https://creativecommons.org/licenses/by/4.0/